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ABSTRACT 


Snow  accumulation  is  a  significant  factor  for  hydrological  planning,  flood  prediction, 
traffic  ability,  avalanche  control,  and  numerical  weather/climatological  modeling. 
Current  snow  depth  methods  fall  short  of  requirements.  This  research  explores  a  new 
approach  for  determining  snow  depth  using  airborne  interferometric  synthetic  aperture 
radar  (InSAR).  Digital  elevation  models  (DEM)  are  produced  for  Snow  Off  and  Snow 
On  cases  and  differenced  to  determine  elevation  change  from  accumulated  snow. 
Interfero grams  are  produced  using  Multi-pass  Single  Look  Complex  airborne  Ku-band 
SAR.  Two  approaches  were  attempted.  The  first  is  a  classical  method  similar  to 
spaceborne  InSAR  and  relies  on  detennining  the  baseline  of  the  interferometric  pair.  The 
second  used  a  perturbation  method  that  isolates  and  compares  high  frequency  terrain 
phase  to  elevation  to  generate  a  DEM.  Manual  snow  depth  measurements  were  taken  to 
verify  the  results.  The  first  method  failed  to  obtain  a  valid  baseline  and  therefore  failed. 
The  second  method  resulted  in  representative  DEMs  and  average  snow  depth  errors  of  - 
8cm,  95cm,  -49cm,  176cm,  87cm,  and  42cm  for  six  SAR  pairs  respectively. 
Furthennore,  Ku-band  appeared  to  be  a  high  enough  frequency  to  avoid  significant 
penetration  of  the  snow.  Results  show  that  this  technique  has  promise  but  still  requires 
more  research  to  refine  its  accuracy. 
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I.  INTRODUCTION 


The  ability  to  measure  snow  depth  over  a  large  area  is  a  capability  that  has  eluded 
the  operational  and  scientific  communities  for  many  years.  This  research  is  a  first  step 
towards  developing  methods  for  detennining  snow  depth  utilizing  synthetic  aperture 
radar  (SAR)  technology.  Because  the  SAR  hardware  is  already  fielded,  successful 
measurement  of  snow  depth  using  remote  sensing  approaches  will  potentially  lead  to 
operational  capabilities  beyond  what  is  currently  available  using  a  variety  of  other 
methods. 

There  are  many  reasons  for  both  military  and  civilian  agencies  to  have  the  ability 
to  quantify  snow  depth.  The  U.S.  military’s  need  to  estimate  the  snow  depth  and  its 
effects  every  spring  in  Afghanistan  in  particular  is  the  motivation  for  pursuing  this 
research  to  develop  a  new  and  operationally  relevant  method  for  detennining  snow  depth. 

While  there  is  a  need  for  the  military  to  have  this  capability  in  Afghanistan,  there 
are  a  multitude  of  other  reasons  that  justify  the  pursuit  of  developing  an  improved 
technique  to  address  this  issue.  The  amount  of  snow  on  the  ground  before  the  spring  melt 
plays  a  critical  role  in  several  areas.  Water  management  is  the  first  and  most  critical  of 
these.  Urban  and  agricultural  communities  can  be  directly  affected  by  the  amount  of 
stored  fresh  water  contained  in  the  seasonal  snow  fall.  This  is  especially  true  in  dry 
climates  that  rely  on  high  altitude  snow  melt  to  replenish  fresh  water  reservoirs.  Central 
and  southern  California  are  prime  examples  of  this.  In  addition  to  water  management,  the 
spring  snow  melt  can  be  a  catalyst  for  swollen  rivers  and  flooding.  An  accurate  estimate 
of  the  snow  volume  in  a  particular  watershed  can  go  a  long  way  towards  forecasting  the 
spring  flood  potential.  This  information  can  be  critical  to  local  authorities  for  risk 
mitigation. 

Trafficability,  both  on  and  off  road,  is  another  area  impacted  by  the  amount  of 
snow  on  the  ground.  Snow  depth  can  detennine  where  and  what  kind  of  operations  can 
be  carried  out  and  if  there  is  a  need  for  special  equipment  to  conduct  these  operations. 
Additionally,  this  information  can  be  used  to  predict  an  adversarial  force’s  movement  and 
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capability  to  carry  out  their  operations.  Trafficability  is  not  restricted  to  just  ground 
movement.  The  ability  to  perform  insertion  or  extraction  missions  by  helicopter  and 
airborne  operations  are  additionally  impacted  by  snow  cover.  For  example,  the  presence 
and  amount  of  snow  can  cause  potential  white-out  conditions,  making  it  impossible  to 
land  rotary  wing  aircraft  for  either  standard  operations  or  medical  evacuation  (medivac) 
emergencies. 

Avalanche  risk  assessment  is  yet  another  area  that  can  benefit  from  improved 
snow  depth  estimates.  One  of  the  many  factors  in  determining  the  potential  for 
avalanches  is  the  amount  of  snow  on  a  particular  slope.  An  estimation  of  the  snow 
volume  along  with  an  estimated  density  can  be  used  in  conjunction,  with  other  factors 
such  as  slope  angle,  seasonal  history,  and  weather  conditions  to  forecast  snow  pack 
stability  (American  Avalanche  Association  2004).  This  generally  is  not  as  big  of  a 
problem  in  areas  that  have  an  active  avalanche  control  system  such  as  those  found  in 
many  recreational  ski  areas.  In  backcountry  areas  and  less  developed  areas  around  the 
world,  however,  snow  depth  can  be  a  valuable  source  of  intelligence  for  both  military 
operations  and  humanitarian  risk  assessment. 

Snow  depth  additionally  plays  a  role  when  it  comes  to  numerical  weather  and 
climate  prediction.  The  high  albedo  of  snow  greatly  impacts  the  local  and  global  scale 
radiational  balance.  The  snow  decreases  the  absorption  of  solar  energy  received  at  the 
earth’s  surface,  acting  as  a  cooling  mechanism  (Groisman  et  al.  1994).  Accurately 
predicting  when  and  where  the  snow  cover  melts  to  expose  bare  ground  benefits 
numerical  weather  and  climatological  model  estimates  of  longwave  and  shortwave 
energy  exchange.  Having  a  more  accurate  estimate  of  the  incoming  energy  will, 
theoretically,  improve  forecasts. 

Current  methods  to  estimate  snow  depth  fall  well  short  of  our  needs.  While  many 
observation  sites  record  the  amount  of  snowfall,  few  record  the  actual  snow  depth.  Those 
that  do,  do  not  come  close  to  providing  adequate  coverage.  Furthermore,  the  few  reports 
of  snow  depths  that  are  available  are  not  necessarily  representative  of  the  local  area  due 
to  spatial  variability. 
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Considering  the  relevance  of  this  infonnation,  the  objective  of  this  research  in  the 
Snow  Depth  Airborne  Radar  (SNODAR)  project  is  to  explore  the  viability  of  a  new 
method  using  remote  sensing,  and  more  specifically  interferometric  synthetic  aperture 
radar  (InSAR)  in  the  Ku-band  (12-18Ghz  frequency  range)  from  an  airborne  platform,  to 
determine  snow  depth  in  regions  that  are  not  easily  accessible.  This  new  method  should 
be  able  to  provide  both  higher  spatial  resolution  and  accuracy  for  point  snow  depth 
measurements  and  area  snow  volume  estimates  than  the  current  methods  available.  It  is 
anticipated  that  existing  SAR  hardware  could  be  used  to  develop  the  capability  to 
detennine  both  snow  depth  and  snow  volume  where  manual  measurements  are  not 
possible. 
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II.  BACKGROUND 


This  section  covers  four  different  existing  approaches  that  have  been  taken  to 
address  the  snow  depth  issue  along  with  their  strengths  and  weaknesses.  The  primary 
methods  that  have  been  used  to-date  include  Air  Force  Weather  Agency’s  Snow  Depth 
and  Sea  Ice  Analysis  (SNODEP)  model,  the  use  of  NASA’s  SIR-C/X-SAR  missions,  the 
use  of  ground  penetrating  radar,  and  the  use  of  LiDAR.  The  final  part  of  this  section 
focuses  on  the  method  investigated  by  this  research,  summarizing  the  theory  used  to 
derive  digital  elevation  models  (DEM)  through  the  use  of  airborne  or  spacebome 
synthetic  aperture  radar  (SAR),  and  describing  the  science  and  techniques  used  to 
determine  snow  depth  with  SAR  systems. 

A.  SNODEP 

The  Air  Force  Weather  Agency’s  SNODEP  model  is  the  primary  tool  used  today 
to  provide  military  operational  users  with  snow  depth  information.  Snow  depth  estimates 
are  modeled  by  using  a  combination  of  passive  microwave  imagery  from  the  Special 
Sensor  Microwave/Imager  Sounder  (SSM/IS)  and  surface  observations  to  include 
synoptic,  meteorological  reporting  observations  (METAR)  and  Airways  and  snow  depth 
climatology  (Air  Force  Weather  Agency  [AFWA]  fact  sheet  and  AFWA  Algorithm 
Description  Document  [ADD]  2012). 

SNODEP  takes  a  multiple  source  approach  to  detennine  snow  depth.  The  first 
step  makes  an  initial  guess  based  on  the  previous  model  run,  similar  to  the  approach  used 
in  many  numerical  weather  prediction  models  to  establish  an  initial  background  field. 
Once  the  background  field  is  established,  the  model  incorporates  any  available  surface 
depth  observations.  It  uses  an  inverse  linear  weighting  scheme  to  interpolate  the  data  to 
the  closest  grid  point.  Then,  in  regions  without  surface  reports,  SSMI/S  algorithms  are 
used  to  detect  snow.  If  no  snow  was  previously  detected,  a  value  of  0.1m  of  snow  depth 
is  automatically  assigned.  If  snow  is  detected  where  snow  was  previously  detected,  the 
snow  depth  estimate  is  trended  toward  climatology.  If  no  snow  is  detected,  the  estimate 
for  the  area  remains  snow  free  (AFWA  fact  sheet  and  AFWA  ADD  2012). 
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The  main  strength  of  SNODEP  is  its  ability  to  provide  a  global  view  of  snow 
coverage.  It  does,  however,  have  several  weaknesses.  Due  to  the  inherent  resolution  of 
the  SSM/IS  satellite;  SNODEP’s  best  resolution  is  25km  (Foster  2011).  This  spatial 
resolution  typically  is  not  adequate  to  provide  the  detail  that  operational  users  require.  Its 
grid  can  also  be  too  large  to  adequately  estimate  the  snow  depth  in  a  smaller  watershed, 
especially  in  complex  terrain  such  as  mountainous  regions.  In  addition,  the  in-situ 
observations  are  extremely  limited.  AFWA  estimates  that  between  100  and  150  snow 
depth  observations  across  the  northern  hemisphere  in  winter  are  used  to  feed  the  model 
(T.  Lewiston  2012,  personal  communication).  This  number  of  observations  does  not 
begin  to  provide  a  good  estimate  for  the  generation  of  a  hemispheric-scale  model.  Also, 
the  observations  that  we  do  have  tend  to  be  concentrated  in  more  developed  countries  like 
the  U.S.  This  makes  it  unlikely  that  these  observations  will  have  an  effect  on  snow  depth 
estimates  for  locations  where  military  operations  are  taking  place.  Spatial  variability  is 
also  a  big  concern.  Wind,  complex  terrain,  and  variable  solar  radiation  are  the  main 
causes  of  this  variability  (Schweizer  and  Kronholm  2006).  Because  of  this  high  spatial 
variability  in  snow  depth,  the  observations  that  are  available  may  not  be  representative  of 
the  region  that  SNODEP  is  trying  to  describe.  To  make  up  for  this  poor  coverage  of  in- 
situ  observations  the  SSMIS  passive  microwave  satellite  is  used  to  detennine  the  snow 
depth  everywhere  else.  SSMIS  does  this  by  using  a  correlation  coefficient  between  the 
microwave  brightness  temperature  and  snow  depth.  This  coefficient  assumes  snow 
crystal  grain  size,  and  that  the  snow  is  dry  or  refrozen.  Failure  of  either  of  these 
assumptions  can  negatively  affect  the  accuracy  of  the  model.  Furthennore,  snow  depth 
estimates  from  the  SSMIS  are  limited  to  depths  of  40cm  or  less.  The  snow  depth 
algorithm  becomes  unreliable  when  the  snow  depth  exceeds  40cm  (Northrop  Grumman 
2010).  Additionally,  at  first  indication  of  snow  on  the  ground,  the  model  automatically 
assigns  0.  lm  of  snow  depth  (AFWA  ADD  2012).  This  may  or  may  not  be  representative 
of  the  initial  snow  depth.  When  snow  was  previously  detected  by  the  algorithm  it  trends 
the  current  value  towards  climatology  (AFWA  ADD  2012).  There  are  two  problems  with 
this.  The  first  is  that  the  current  snow  depth  may  or  may  not  be  representative  of  the 
snow  on  the  ground  at  the  modeled  time.  Changing  the  value  is  premature  if  the  depth  is 
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not  accurate  to  begin  with.  The  second  problem  has  to  do  with  the  climatological  data. 
As  mentioned  above,  there  are  relatively  few  stations  that  actually  record  snow  depth. 
Many  stations  record  snowfall,  which  should  not  be  confused  with  snow  depth  on  the 
ground.  Mechanisms  such  as  settling,  melting,  sublimation,  and  movement  of  snow  by 
wind  make  the  snowfall  measurements  a  poor  estimate  of  snow  depth.  With  so  few 
stations  recording  snow  depth,  the  accuracy  of  the  climatology  is  in  question,  especially 
for  areas  that  do  not,  or  never  have  had,  in-situ  measurements.  This  potentially  can  have  a 
negative  impact  on  the  accuracy  of  SNODEP  estimates. 

B.  SIR-C/X-SAR 

1.  SIR-C/X-SAR  Description 

The  Spacebome  Imaging  Radar-C/X-band  Synthetic  Aperture  Radar  (SIR-C/X- 
SAR)  mission  was  a  joint  National  Aeronautical  and  Space  Administration  (NASA), 
Gennan  Space  Agency  (Deutsche  Agentur  fur  Raumfahrtangelegenheiten,  DARA),  and 
Italian  Space  Agency  (Agenzia  Spaziale  Italiana,  ASI)  project  designed  as  the  next  step 
in  spaceborne  imaging  radar  research.  SIR-C/X-SAR  flew  two  missions  in  1994  imaging 
57.6  million  square  miles,  or  approximately  14  percent  of  the  Earth’s  surface  (Stofen  et 
al.  1995,  and  JPL,  2012).  The  space  shuttle  launched  with  three  different  synthetic 
aperture  radar  (SAR)  antennas.  These  included  L-band  (23.5cm  wavelength),  C-band 
(5.8cm  wavelength),  and  X-band  (3cm  wavelength)  antennas.  The  L  and  C  bands  were 
also  constructed  in  a  way  that  could  measure  both  horizontal  and  vertical  polarizations. 
The  use  of  the  three  different  bands  allowed  collection  of  information  about  the  Earth’s 
surface  at  multiple  scales,  which  had  never  been  possible  before  with  only  single  band 
SAR  systems. 


2.  Snow  Analysis  Using  SIR-C/X-SAR 

Snow  characteristics  greatly  affect  backscattering  of  radar  emissions.  The  extent 

of  this  effect  is  dependent  on  three  sets  of  parameters  (Shi  and  Dozier  1996): 

•  Sensor  parameters  to  include  frequency/wavelength,  polarization,  and 
viewing  angle 
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•  Snow  pack  parameters  to  include  snow  density,  depth,  particle  size,  size 
variation,  liquid  water  content  (stickiness),  and  stratification 

•  Ground  parameters  to  include  dielectric  and  roughness  parameters 

This  difference  in  backscattering  properties  by  different  radar  wavelengths  on  the 
snow  pack  can  be  leveraged  to  detennine  the  physical  characteristics  of  the  snow  pack 
and  the  underlying  ground.  All  three  of  the  SIR-C/X-SAR  wavelengths  are  assumed  to 
penetrate  into  the  snowpack.  Based  on  electro-magnetic  scattering  theory,  for  a  given 
material,  there  is  a  direct  relationship  between  the  wavelength  and  the  depth  of 
penetration  (Richards  2009).  With  that  in  mind,  there  should  be  an  increase  in 
backscattering  moving  from  the  L-band  radar  down  to  the  X-band  radar.  This  fact  was 
used  by  Shi  and  Dozier,  2000,  to  retrieve  snowpack  properties.  They  first  used  polarized 
data  from  the  L-band  radar  to  determine  snowpack  density.  L-band  proved  to  be  a  long 
enough  wavelength  that  the  backscatter  from  the  snowpack  was  negligible.  The  entire 
radar  return  therefore  came  from  the  ground  below  the  snow  pack.  Despite  the  lack  of 
backscatter  from  the  snow  they  were  able  to  capitalize  on  the  fact  that  the  snow  pack 
caused  a  shift  in  refraction  in  the  incidence  angle  of  the  radar  pulse.  The  extent  of  the 
refraction  was  dependent  on  the  density  of  the  snow  pack.  Furthermore,  there  was  a 
difference  in  both  the  magnitude  and  relation  between  the  VV  and  HH  polarizations.  By 
modeling  this  interaction,  they  were  able  to  derive  the  snowpack’s  density. 

Due  to  the  large  variability  in  density  in  snowpacks,  the  density  alone  is  not 
enough  to  estimate  other  characteristics  of  the  snow  pack  such  as  snow  depth  or  snow 
water  equivalent  (SWE).  To  do  this  Shi  and  Dozer  used  data  from  both  the  C-band  and 
X-band  radars.  Both  C-band  and  X-band  radar  pulses  have  different  volume  scattering 
properties  as  stated  above.  This  fact  can  be  used  to  model  the  particle  size  and  expected 
magnitude  of  the  scattering.  Both  bands  are  assumed  to  penetrate  to  the  ground  in 
addition  to  the  volume  scattering,  which  adds  an  additional  component  to  the  overall 
return.  This  can  be  accounted  for,  however,  by  using  the  ground  roughness  and  dielectric 
properties  determined  from  the  L-band  radar. 

This  approach  of  using  a  combination  of  all  three  SAR  bands  showed  very 
positive  results  and  has  stood  up  well  to  ground  validation.  While  this  technique  has 
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shown  great  potential,  we  do  not  currently  have  any  persistent  spacebome  or  airborne 
sensors  with  the  appropriate  configuration  to  take  advantage  of  this  technique. 


C.  GROUND  PENETRATING  RADAR 

There  have  also  been  attempts  to  use  ground  penetrating  radar  (GPR)  to  address 
the  issue  of  detennining  snow  depth  and  other  snowpack  characteristics.  Frequency 
modulated  continuous  wave  (FMCW)  radar  has  proven  to  be  the  most  successful  of  the 
GPRs  for  snow  study.  The  FMCW  radar  works  similarly  to  a  standard  radar  system  in 
that  it  times  the  pulse  to  detennine  range.  It  however  uses  a  broad  band  width  that  results 
in  a  greater  theoretical  vertical  resolution  as  compared  to  a  standard  GPR  (Yankielun  et 
al.  2004).  This  greater  vertical  resolution  is  quite  important  if  you  want  to  detennine 
snow  pack  stratigraphy,  which  can  be  particularly  important  for  avalanche  prediction. 
Ground  penetrating  radars  are  typically  deployed  for  snow  pack  analysis  either  by  hand 
or  by  towing  them  behind  a  snowmobile.  Recently  they  have  been  deployed  using  low 
flying  helicopters  with  some  success  (Marshal  et  al.  2008). 

Overall,  the  use  of  these  FMCW  radars  has  been  quite  successful  at  determining 
snowpack  characteristics;  in  particular  those  characteristics  that  concern  avalanche 
experts  in  focused  areas.  They  are  not,  however,  suited  for  covering  larger  areas. 
Deploying  them  on  the  ground,  whether  by  hand  or  being  towed  behind  a  snowmobile  or 
snowcat,  does  not  provide  nearly  the  spatial  coverage  provided  by  airborne  systems. 
Ground  deployment  is  also  restricted  by  complex  terrain.  The  use  of  the  GPR  by 
helicopter  also  has  drawbacks.  The  systems  used  to-date  have  a  fairly  broad  footprint. 
That  means  that  as  the  platfonn  that  the  GPR  is  attached  to  increases  in  altitude,  the  area 
covered  by  the  footprint  also  increases  dramatically.  Everything  in  the  footprint  is  treated 
as  a  single  return  per  pulse.  The  more  the  terrain  varies  within  the  footprint,  the  less 
reliable  the  measurements.  Work  done  by  Marshal  et  al.  (2008)  has  shown  that  altitudes 
greater  than  100ft  above  the  ground  make  the  data  unreliable.  Performance  can  be  worse 
in  areas  where  there  are  steep  slopes.  There  are  plans  to  try  to  use  a  FMCW  GPR  with  a 
narrower  beam  to  address  this  issue.  At  the  time  of  this  writing,  however,  there  do  not 
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appear  to  be  publications  indicating  that  this  has  been  done  yet.  With  such  restrictions, 
operational  airborne  collections  in  complex  terrain  are  not  currently  possible. 

D.  LIDAR 

Light  Detection  and  Ranging  (LiDAR)  is  another  method  that  has  been  explored 
to  estimate  snow  depth.  LiDAR  is  based  on  measuring  the  time  required  for  a  pulse  of 
light  to  travel  to  a  target  and  then  return  to  determine  range  (Hodgson  et  al.  2005).  This 
can  be  used  to  build  either  2-dimensional  or  3-dimensional  scenes.  To  detennine  snow 
depth,  the  scene  is  imaged  with  and  without  snow  and  then  differenced,  resulting  in  a 
snow  volume  and  snow  depth  estimate  at  each  specific  point.  The  use  of  LiDAR  has  a  lot 
of  advantages.  The  first  is  that  it  can  be  used  to  cover  large  areas  in  an  unobtrusive 
manner.  It  is  also  highly  accurate,  with  accuracies  down  to  the  millimeter  level  in  some 
cases  (Osterhuber  et  al.  2008). 

LiDAR  has  been  deployed  two  different  ways  to  determine  snow  depth.  The  most 
accurate  way  is  to  deploy  the  LiDAR  system  on  the  ground.  Osterhuber  et  al.  (2008) 
used  a  ground  based  unit  that  could  either  be  placed  on  the  ground  or  fixed  to  a 
surveyor’s  tripod.  In  a  snow  pack  with  an  average  depth  of  just  over  two  meters  the 
LiDAR  averaged  a  mean  difference  between  manual  and  LiDAR  measurements  of  5.7 
cm.  While  the  use  of  the  ground-based  system  has  potential,  it  also  has  some  drawbacks. 
The  system  currently  being  used  has  a  range  limited  to  1000  meters.  Also,  to  generate  a 
3-D  image,  either  multiple  sensors  are  required  or  the  LiDAR  has  to  be  moved  to 
different  scanning  locations.  Furthermore,  LiDAR  becomes  ineffective  with  any 
obscuring  weather  phenomena  such  as  clouds,  fog,  or  precipitation.  This  system  may 
prove  to  be  a  great  way  of  measuring  snow  depth  at  fixed  locations  but  is  not  a  good 
option  for  large  regions  of  land  or  remote  areas  where  a  ground-based  unit  has  not  been 
placed. 

The  second  way  to  deploy  LiDAR  is  to  operate  the  system  from  either  an  aircraft 
or  a  spaceborne  platform.  Airborne  LiDARs,  also  known  as  laser  altimetry,  are  much 
better  suited  to  cover  large  regions  or  remote  areas  than  the  fixed  based  systems 
(Hodgson  et  al.  2005).  Airborne  LiDAR  depends  on  knowing  the  speed  of  light,  the 
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location  of  the  laser  emitter,  and  being  able  to  time  the  laser  pulse  transmission  to 
reception  time.  These  data,  like  the  ground  based  systems,  can  be  used  to  generate  a  3-D 
image  or  terrain  model  with  a  resolution  at  sub-meter  level  (Hopkinson  et  al.  2004).  This 
has  the  same  restriction  as  the  ground  based  system  in  the  fact  that  the  laser  path  has  to  be 
free  of  visual  obscurations.  Accuracy  also  depends  on  the  ability  to  position  the  aircraft 
to  a  high  degree  of  x,  y,  z  accuracy,  which  can  potentially  be  problematic.  Furthermore, 
there  are  a  limited  number  of  platforms  that  are  currently  equipped  to  perform  this  task. 

E.  INSAR  AND  TOMOGRAPHIC  SAR 

The  approaches  investigated  in  this  research  are  based  on  radar  (radio  detection 
and  ranging)  characteristics,  principals  and  theory,  utilizing  Interferometric  Synthetic 
Aperture  Radar  (InSAR).  Discussion  of  this  topic  first  requires  a  review  of  Synthetic 
Aperture  Radar  (SAR)  basics.  Radar  uses  radiation  emitted  from  an  antenna  in  the 
microwave  region  of  the  electromagnetic  (EM)  spectrum.  The  wavelengths  most 
commonly  used  are  on  the  order  of  1.5cm  tolm,  or  approximately  20GHz  to  300MHz 
(Richards  2009)(Figures  1  and  2).  This  frequency  range  is  broken  down  into  bands  with 
L-  (1. 0-2.0  GHz),  C-  (4. 0-8.0  GHz),  and  X-bands  (8.0-12.5  GHz)  as  the  most  commonly 
used  for  remote  sensing  as  seen  in  Figure  2.  This  emitted  energy  travels  to  a  target  and  is 
then  reflected  back  to  the  original,  or  in  some  cases  an  alternate  antenna.  The  time  it  takes 
this  radiation  to  travel  the  distance  to  and  then  back  from  the  target  is  measured.  Using 
the  speed  of  electromagnetic  propagation,  this  allows  an  estimate  of  the  range  to  the 
target  (Carrara  et  al.  1995).  The  properties  of  the  microwave  portion  of  the  EM  spectrum 
in  particular  make  it  a  good  option  for  transmission  through  the  atmosphere.  Microwave 
frequencies  are  able  to  pass  freely  through  the  atmosphere  and  through  visible 
obscurations  such  as  clouds  and  precipitation,  allowing  them  to  detect  objects  where 
optical  wavelengths  would  be  either  absorbed  or  scattered.  The  atmospheric 
transmittance  per  frequency  is  annotated  in  Figure  1.  Observe  that  the  lower  end  of  the 
of  the  microwave  region  used  for  remote  sensing,  or  about  20GHz,  has  strong 
transmittance  through  the  atmosphere  and  rapidly  increases  to  near  100  percent  as  the 
frequency  decreases. 
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Figure  1.  The  electromagnetic  spectrum  and  the  indicative  transmittance  of  the  atmosphere 

on  a  path  between  space  and  the  earth  (Richards  2009). 
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Figure  2.  The  microwave  portion  of  the  electromagnetic  spectrum  is  broken  down  into 
bands  with  L-  (1. 0-2.0  GHz),  C-  (4. 0-8.0  GHz),  and  X-bands  (8.0-12.5  GHz)  as 
those  most  commonly  used  for  remote  sensing  (Richards  2009). 


SAR  takes  the  use  of  radar  one  step  further.  Typically  SAR  systems  are 
configured  in  a  side-looking  manner.  Real  aperture  radar  (RAR)  has  limitations  in 
resolution  based  on  the  beam-width  “p”  and  the  distance  to  the  target.  The  beamwidth  is 
dependent  on  the  wavelength  “A,”  and  the  physical  length  of  the  antenna  “la”  and  is  given 
by  equation  (1)  (Jakowatz,  et  al.  1996) 


P  = 


X 


(1) 


When  the  radar  pulse  travels  away  from  the  radar  it  spreads  out.  This  beam 
spreading  can  be  seen  in  Figure  3.  The  radar  antenna  will  receive  a  return  from  all 
objects  within  the  footprint  that  are  at  an  equal  distance  from  the  source.  Because  the  arc 
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of  the  radar  increases  with  increased  range,  the  azimuthal  distance  required  to  distinguish 
between  targets  will  also  increase  with  range.  In  other  words,  the  ability  to  resolve 
objects  of  a  certain  size  in  the  azimuth  direction  of  a  side-looking  radar  will  decrease  as 
distance  from  the  antenna  increases.  This  can  be  a  real  problem  when  using  radar  for 
imaging  rather  than  for  detection.  The  azimuthal  resolution,  “ra”,  on  the  ground  with  a 
range  Ro  and  beam- width  P  is  given  by  (2)  (Jakowatz,  1996  and  Richards,  2009). 

ra  =  KP 
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Figure  3.  Imaging  scenario  depicting  the  cross-range  resolvability  issue.  Targets  lying  on 
the  same  constant-range  line  but  separated  in  cross  range  cannot  be  distinguished 
in  the  return  of  a  single  pulse  when  a  broad  beam  is  used  to  illuminate  the  entire 
ground  patch,  where  L  represents  the  patch  radius  and  p(y)  represents  range  on 
the  ground  in  the  direction  perpendicular  to  the  flight  path  (Jakowatz,  1996). 

As  mentioned  in  the  preceding  paragraph,  the  beam  width  is  dependent  on  the  size 
of  the  antenna.  The  larger  the  antenna,  the  smaller  the  beam  width,  and  the  greater  the 
resolution.  Unfortunately  there  are  limitations  to  the  physical  size  of  the  antennas  that 
can  be  used,  especially  for  airborne  and  spaceborne  radars.  SAR  solves  this  problem  by 
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using  the  forward  motion  of  the  platform  to  which  the  radar  is  attached  to  synthesize  an 
apparently  longer  antenna,  as  seen  in  Figure  4  (Richards  2009).  SAR  irradiates  a  target 
multiple  times  as  it  passes  that  target.  Through  processing  of  the  multiple  images,  a 
larger  antenna  or  narrower  beam  is  synthesized.  The  SAR  azimuthal  resolution  ra  is 
given  by  (3)  where  physical  antenna  length  is  “la” 


r  =  — 
2 


(3) 


Notice  that,  unlike  real  aperture  radar,  there  is  no  dependence  on  the  absolute 
range  from  the  antenna  to  the  target,  or  the  wavelength  used.  This  is  what  allows  SAR  to 
have  improved  resolution  at  the  ground  compared  to  RAR  (Richards,  2009). 
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Figure  4.  The  concept  of  using  the  platfonn  motion  to  synthesize  an  effectively  long 

antenna;  the  footprint  of  the  real  antenna  on  the  ground  is  shown  as  rectangular 
for  simplicity,  where  “La”  represents  the  length  of  the  synthetic  antenna,  “la” 
represents  the  physical  length  of  the  antenna,  “Ro  ”  represents  the  range,  “A” 
represents  the  wave  length  and  “0a”  is  the  beamwidth  (Richards  2009). 


As  a  side  note,  there  is  no  significant  distinction  between  SAR  and  RAR  with 


respect  to  slant  range  resolution.  Ground  range  (the  range  directly  below  the  antenna  to 
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the  target)  resolution  is  dependent  on  the  slant  range  resolution.  Both  the  slant  range  and 
ground  range  resolutions  are  dependent  on  the  width  of  the  radar  pulse  “x”  and  the 
incidence  angle  “0”  (Figure  5).  The  slant  range  resolution  “rr”  can  be  given  by 

CT 


r  - 


2 


The  speed  of  light  “c”  is  treated  as  a  constant, 
ground  range  resolution  “rg”  can  be  given  as 


(4) 

Correcting  for  the  incident  angle,  the 


CT 

s  2  sin  9 


(5) 


the  angle  that  the  incoming 
beam  makes  with  the 
normal  to  the  surface  is 


separation  of  targets  in 
the  ground  range  direction 


Figure  5.  Geometry  for  computing  range  resolutions.  “rg  ”  represents  the  difference  in 
ground  range  with  respect  to  the  radar  antenna,  “rr”  is  the  slant  range  resolution, 
“Ar”  the  difference  in  the  slant  range  between  targets  A  and  B,  and  “0”  is  the 
incidence  angle  between  the  incoming  beam  and  the  nonnal  to  the  surface. 

(Richards  2009) 


Interferometric  SAR  (InSAR)  is  a  quantitative  step  beyond  SAR  imaging.  InSAR 


capitalizes  on  the  capability  to  measure  the  phase  angle  of  the  SAR  return.  The 
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transmitted  phase  is  known  and  the  return  phase  can  be  measured.  This  allows 
detennination  of  relative  distances  from  the  sensor  to  the  ground,  and  thus  the 
detennination  of  topography  or  topographic  displacement.  InSAR  has  been  used  for  the 
last  two  decades  for  topographical  mapping  or  DEM  generation  (Richards  2009). 

The  aforementioned  difference  between  the  transmitted  phase  and  the  received 
phase,  (j>T,  is  given  by  (6)  (Richards  2009) 

,  4  ttR 


where  “R”  represents  the  one  way  path  length  and  “A.”  is  the  operating  wavelength  of  the 
radar.  Observe  that  the  phase  difference,  4>t.  will  change  between  two  adjacent  pixels  as 
the  total  path  distance  changes.  In  other  words  the  phase  angle  change  is  proportional  to 
the  difference  in  path  lengths  between  targets.  Therefore,  since  we  can  measure  phase  we 
should  be  able  to  directly  measure  the  difference  in  the  path  length  between  targets  to  a 
very  high  degree  of  accuracy.  There  is  however  a  potential  ambiguity.  Any  change  in 
elevation  of  an  object  that  stands  above  or  below  the  plane  being  imaged  will  cause  a 
coincident  return  of  the  same  phase.  This  can  be  seen  in  Figure  6  (a)  between  points  “a” 
and  “c”.  This  ambiguity  prevents  us  from  being  able  to  resolve  the  difference  between 
“a”  and  “c”. 
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Figure  6.  (a)  Even  though  the  topographic  variation  between  a  and  b  is  resolvable,  there  is 

ambiguity  between  a  and  c  (b)  resolving  the  a  and  c  ambiguity  by  changing  the 
radar  position,  but  causing  ambiguity  at  a  and  d  (Richards  2009). 


This  ambiguity  can  be  overcome  by  changing  the  position  of  the  radar  as  seen  in 
(b).  The  change  in  the  radar’s  position  however  causes  other  potential  ambiguities,  as 
can  be  seen  between  points  “a”  and  “d”.  Considering  both  view  angles  at  the  same  time, 
the  ambiguities  can  be  resolved.  This  is  the  basis  of  interferometry  (Richards,  2009), 
(Figure  7). 
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Figure  7.  Geometry  for  single  baseline  SAR  interferometry,  in  which  the  look  and 

incidence  angles  are  assumed  to  be  the  same.  “Ri”  and  “R2”  are  the  respective 
ranges  from  antennas  1  and  2.  “B”  represents  the  baseline  between  the  two 
antenna  locations.  “BA  represent  the  orthogonal  baseline  between  the  two  radar 
beam  paths.  “0”  and  “50”  represent  the  incidence  angle  and  the  change  in 
incidence  angle  respectively  (Richards  2009). 

Figure  7  and  accompanying  discussion  illustrates  how  these  two  view  angles  can 
be  used  and  introduces  the  concept  of  the  orthogonal  baseline.  Points  1  and  2  indicate  the 
positions  of  the  platform  where  the  two  datasets  are  acquired.  Line  “B”  that  connects 
these  two  positions  indicates  the  baseline.  A  line  is  drawn  perpendicular  from  the  first 
slant  range  to  position  2  and  is  referred  to  as  the  orthogonal  baseline,  “B  ”.  There  are 
limitations  on  the  length  of  the  orthogonal  baseline  known  as  the  critical  baseline.  This 
is  important  because  it  allows  the  measurement  of  the  difference  in  the  length  of  the  slant 
range  in  terms  of  the  phase  or  as  a  percentage  of  the  wavelength.  This  change  in  phase 
between  the  different  view  angles  of  the  same  target  is  how  InSAR  can  be  used  to 
generate  DEMs.  This  process  is  also  referred  to  as  tomographic  SAR. 

Using  Figure  7,  the  difference  in  the  path  lengths  “Ri”  and  “R2”  in  terms  of  the 
phase  and  a  given  baseline  and  incidence  angle  of  “B”  and  “0”  respectively  as  shown  by 
(Richards  2009)  can  be  derived. 
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i?,  =  R2  cos  SO  +  B  sin  6 


50  is  assumed  to  be  approximately  0  using  the  plane  wave  approximation.  The  plane 
wave  approximation  considers  the  change  in  the  incidence  angle  to  approximate  0  when 
the  target  is  infinitely  far  away  when  compared  to  the  length  of  orthogonal  baseline.  This 


results  in 


Therefore 


Rl  =  R2  +  B  sin  0 


A R  -  Rx—  R2  —  B  sin  6 


The  difference  in  phase  angle  “A(|)”  associated  with  the  change  in  path  length  “AR” 

between  the  two  passes  can  then  be  given  as 

.  4  nB  sin  6 

AS  = - 

A  (10) 

This  difference  in  phase  angle  is  referred  to  as  interferometric  phase  angle  Ac)).  A«j)  can 
be  obtained  directly  by  simply  imaging  an  area  twice  and  taking  the  difference  of  the  two 
recorded  phases. 

The  next  step  is  to  determine  the  relationship  between  the  topographic  height  “h” 
and  the  incidence  angle  in  order  to  get  the  phase  to  height  ratio  (Figure  8). 

From  Figure  8,  if  “H”  is  the  total  height  above  an  assumed  altitude,  and  “Ro”  is  the  range 
to  the  target,  observe  that 

h  =  H  -  R0  cos  6 

Taking  the  partial  derivative  of  the  topographic  height  with  respect  to  the  incidence  angle 
results  in 


R0  sin  6 


Then  taking  the  partial  derivative  of  the  interferometric  phase  angle  Acj)  with  respect  to 
the  incidence  angle  also  results  in 

d{A(j>)  _  4/zScos# 

dO  A  (13) 
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Figure  8.  Determining  the  relationship  between  topographic  height  “h”  and  incidence  angle 
“0”  with  a  platform  altitude  of  “H”  and  range  to  the  target  of  “Ro”  (Richards 

2009). 


Combining  equations  (12)  and  (13)  results  in 

d{A(f))  d(A(f>)  dO  AnB  cos  6 


(14) 


dh  dO  d(h)  /li?osin0 
We  now  have  an  expression  for  the  change  in  interferometric  phase  with  respect  to  the 
change  in  topographic  height.  Taking  it  one  step  further  to  make  it  more  user  friendly 
results  in 

d(A</>)  _  AkBl  _  AkBlcos6 
jo  a  in~r  _  a 

(15) 


dh  ARq  sin  0  A(H  -  h )  sin  0 

So  as  long  as  the  incidence  angle  is  known,  the  elevation  above  some  known  reference 


height  (H-h)  and  the  orthogonal  baseline,  the  rate  of  change  in  elevation  across  an 
interferometric  phase  diagram  per  change  of  radian  can  be  predicted.  An  interferometric 
phase  factor  <Xif  can  be  defined  as 

dh 


aiF 


d(A</>) 


and  the  height  of  a  specific  pixel  will  be  given  by 
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(16) 


h(x,y)  =  aIFA<f>(x,  y)  +  CONSTANT 


(17) 

Equation  (17)  enables  the  generation  of  a  DEM.  This  ability  to  use  InSAR  to 
generate  a  DEM  is  ultimately  the  basis  of  this  research.  The  use  of  high  resolution  DEMs 
will  be  used  to  explore  a  new  technique  to  solve  for  the  snow  depth. 
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III.  METHODOLOGY 


A.  RESEARCH  APPROACH 

This  research  used  SAR  data  in  an  approach  similar  to  that  taken  for  airborne 
LiDAR  determination  of  DEMs.  Interferometric  (tomographic)  SAR,  InSAR,  was  used 
instead  of  laser  altimetry  to  map  both  the  bare  ground  and  the  snow  covered  ground. 
Snow  cover  effectively  acts  to  change  the  elevation  of  the  surface.  Taking  the  difference 
between  the  two  DEMs,  Snow  On  and  Snow  Off,  allows  for  detennination  of  snow  depth 
and  snow  volume  over  a  specified  area.  InSAR  has  a  couple  of  advantages  over  LiDAR. 
Because  SAR  operates  in  the  microwave  region  of  the  electromagnetic  spectrum,  it  can 
be  used  when  there  are  obscurations  such  as  clouds  and  precipitation.  These  obscurations 
are  common  during  winter  and  can  be  a  limiting  factor  for  the  laser  based  systems.  From 
an  operational  standpoint,  InSAR  has  another  advantage.  There  are  numerous  platforms 
currently  carrying  SAR  for  other  purposes,  most  notably  the  MQ-1  Predator  and  MQ-9 
Reaper  (General  Atomic  Aeronautical  2012).  These  can  potentially  be  adapted  for 
operational  InSAR  snow  depth  detennination. 

Two  different  SAR  approaches  were  explored.  The  first  approach,  refened  to  as 
the  “Baseline  Method”,  attempted  to  calculate  the  baseline  between  two  different  data 
sets’  locations.  The  second  method,  refened  to  as  the  “Best  Fit  Plane  Removed  Method” 
(BFPR),  circumvents  the  need  to  calculate  the  baseline  and  determines  the  difference  in 
the  Snow  On  and  Snow  Off  terrain  by  removing  a  best  fit  plane  (BFP). 

The  proposed  approaches  have  not  previously  been  taken  before  for  snow  depth 
measurement  for  a  couple  of  reasons.  Generally  speaking,  active  airborne  and 
spaceborne  radar  tend  to  penetrate  the  snowpack,  making  it  impossible  to  infer  the  snow- 
air  interface.  This  is  primarily  due  to  the  frequency  (wavelengths)  used.  Most  SAR 
systems  tend  to  use  L-  (1. 0-2.0  GHz),  C-  (4. 0-8.0  GHz),  and  X-band  (8.0-12.5  GHz) 
radars.  All  of  these  bands  significantly  penetrate  the  snow  pack  (Ulaby,  1982).  In  these 
bands,  with  the  exception  of  using  a  radar  with  a  relatively  high  power  output  located  at 
close  range  to  a  target  such  as  that  used  with  ground  penetrating  radar,  there  is  not 


25 


enough  backscatter  from  the  snow  pack  to  identify  the  snow-air  interface.  The  longer 
wavelengths  primarily  show  volume  scattering,  which  varies  depending  on  the  condition 
of  the  snow.  It  has  been  demonstrated  by  Ulaby,  1982,  however,  that  shorter  wavelength 
(higher  frequency)  SAR  systems  such  as  Ku-  and  Ka-band  do  not  penetrate  the  snow 
pack  as  far  as  other  longer  wavelength  systems,  Figure  9.  This  research  takes  advantage 
of  the  fact  that  the  SAR  used  on  the  MQ-1  Predator  and  MQ-9  Reaper  platforms  uses  a 
shorter  wavelength  or  higher  frequency  (Ku-band)  radar.  This  higher  frequency  should 
penetrate  less  into  the  snow  pack,  resulting  in  a  greater  interaction  at  the  snow-air 
interface  and  leading  to  successful  snow  depth  determination. 


SNOW 


Volumetric  Liquid  Water 
Content  mv  (Percent) 

Figure  9.  Penetration  depth  as  a  function  of  frequency.  Ulaby  et  al.  1982  demonstrates  that 
as  frequency  increase,  the  penetration  into  the  snowpack  decreases.  In  contrast, 
lower  frequencies  will  penetrate  further  into  the  snowpack. 
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B.  SAR  DATA  COLLECTION  AND  FIELD  WORK 

Unfortunately,  no  archived  data  sets  were  available  at  the  start  of  this  research 
that  would  allow  exploration  of  the  InSAR  snow  depth  mapping  concept  using  short- 
wavelength  SAR.  Therefore,  in  order  to  explore  this  approach,  a  new  SAR  data  set  was 
collected  and  processed.  Furthennore,  direct  measurements  of  snow  depth  were  made  to 
validate  and  determine  the  accuracy  of  the  results.  The  field  work  and  data  collection 
consisted  of  four  parts;  (1)  aircraft  SAR  data  collection,  (2)  site  selection  and  preparation, 
(3)  site  survey,  and  (4)  validation. 

1.  Aircraft  SAR  Collection 

The  first  portion  of  the  field  work  and  the  ultimate  tool  used  in  this  approach  was 
the  airborne  SAR  collection.  General  Atomics  Aeronautical  (GAA),  the  manufacturer  of 
the  Lynx  Ku-band  radar  (Tsunoda  et  al.  1999),  agreed  to  fly  collection  missions  to 
generate  the  raw  radar  datasets  necessary  for  the  technique.  The  collection  for  this 
research  was  done  using  a  Lynx  II  SAR  mounted  to  a  King  Air  aircraft.  As  stated  earlier 
this  is  the  same  ground  imaging  radar  that  is  used  on  many  military  aircraft.  The  radar 
itself  is  a  Ku-band  radar  that  operates  at  a  15.2-1 8. 2GHz  frequency  with  a  wavelength  of 
1.8cm.  The  maximum  slant  range  of  this  system  is  30km.  There  are  multiple  ground 
resolution  options,  to  include  0.1,  0.3,  1.0,  and  3.0m.  Table  1  indicates  the  SAR  images 
used  in  this  study.  The  radar  can  be  operated  in  different  modes  to  include  spotlight, 
strip-map,  and  coherent  change  detection  (Tsunoda  et  al.  1999).  For  the  purposes  of  this 
research,  multiple  passes  were  made  at  both  the  0.1  and  0.3m  resolutions  using  the 
spotlight  mode.  Both  Snow  On  and  Snow  Off  collections  were  flown  on  3  April  2012 
and  13  July  2012  respectively.  The  actual  flights  were  flown  at  a  flight  level  of  5330m 
and  approximately  10km  to  the  south  of  the  collection  area  with  the  radar  looking  north. 
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SAR  Image 
number 

Date 

Time 

Surface  Condition 

01 

3  April  2012 

18:47  Z,  10:47  L 

snow  covered 

02 

3  April  2012 

18:55  Z,  10:55  L 

snow  covered 

03 

3  April  2012 

19:04  Z,  11:04  L 

snow  covered 

04 

3  April  2012 

19:00  Z,  11:00  L 

snow  covered 

21 

13  July  2012 

17:58  Z,  10:58  L 

bare 

22 

13  July  2012 

18:05  Z,  11:05  L 

bare 

Table  1 .  Four  usable  Snow  On  and  two  useable  Snow  Off  SAR  images  at  0.  lm  resolution 
were  acquired  by  General  Atomics  during  the  radar  collection  phase  of  the 
research.  These  images  enabled  six  Snow  On  and  one  Snow  Off  interferometric 

pairs  to  be  calculated. 

2.  Site  Selection  and  Preparation 

The  second  aspect  of  the  field  work  revolved  around  the  site  preparation.  First  an 
appropriate  site  had  to  be  chosen.  There  were  several  requirements  to  consider  in  the  site 
selection.  They  were:  sufficient  snow  depth,  relatively  obstacle  free,  clear  view  of  sky, 
relatively  flat  i.e.  no  steep  slopes  relative  to  the  resolution,  accessibility,  and  within  the 
General  Atomics  flight  radius  centered  in  San  Diego,  California.  It  must  be  noted  here 
that  the  first  criterion  of  sufficient  snow  depth  played  a  major  role  in  site  selection.  Due 
to  the  low  snow  fall  in  the  2011-2012  winter  and  dissertation  time  restrictions  there  were 
not  many  viable  options.  Mammoth  Mountain,  California  (Figure  10)  ultimately  met  the 
requirements  better  than  any  other  location. 
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Figure  10.  Map  depicting  relationship  between  General  Atomics  home  airfield  and  Snow 
study  site  at  Mammoth  Mountain,  CA  (Google  Earth  2012). 


Mammoth’s  higher  elevation  with  respect  to  other  sites  provided  a  greater  and 
longer  lasting  snowpack.  Ski  lifts  made  many  of  the  higher  elevation  areas  accessible.  It 
was  also  within  the  GAA  flight  radius.  In  addition,  it  also  a  Cold  Regions  Research  and 
Engineering  Laboratory  (CRREL)  snow  study  site  situated  there  and  the  Mammoth 
Mountain  Ski  Patrol  was  both  willing  and  interested  in  participating.  This  ultimately 
made  a  big  difference.  Their  local  knowledge  of  the  mountain  and  surrounding  areas 
were  key  to  identifying  good  locations  in  the  time  allotted.  With  their  help,  eight 
possible  collection  locations  were  identified.  SAR  data  were  ultimately  only  collected 
for  one  site,  “Elysian  Fields,”  shown  in  Figure  11. 
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Figure  11.  Google  Earth  3D  perspective  view  of  Elysian  Fields  collection  area  looking 
northward  towards  Mammoth  Mountain,  CA  (Google  Earth  2012). 


The  generation  of  InSAR  images  and  production  of  DEMs  from  airborne  radar 
data  is  not  a  straightforward  process.  One  of  the  things  that  can  be  done  to  increase  the 
accuracy  is  the  use  of  ground  control  points  (GCP),  also  commonly  referred  to  as  tie 
points.  This  can  be  done  using  natural  features  or  manmade  features  or  equipment. 
Manmade  GCPs,  specifically  comer  reflectors  (CR),  however,  provide  several  orders  of 
magnitude  greater  radar  return  than  that  of  the  surrounding  terrain.  Highly  reflective  CRs 
make  it  much  easier  to  identify  specific  locations  in  each  image.  Their  use  removes 
much  of  the  ambiguity  in  the  DEM  generation  process.  In  addition  to  the  generation  of 
the  InSAR  images  and  DEMs,  the  CRs/GCPs  can  be  used  to  warp  the  Snow  On  and 
Snow  Off  DEMs  for  the  differencing  phase  of  the  processing.  Furthermore,  they  can  be 
used  to  extract  information  that  can  be  used  to  determine  the  baseline  between  the  two 
flight  paths  used  in  the  generation  of  an  interferometric  image.  There  are  several 
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different  configurations  that  can  be  used  in  the  construction  of  corner  reflectors. 
Examples  of  these  are  shown  Figure  12. 


Figure  12.  Backscatter  from  various  shapes  and  corresponding  radar  cross  section  equations 

(Granit  Island  Group,  2013). 


With  this  in  mind,  trihedral  radar  reflectors  were  chosen  to  be  deployed  at 
selected  ground  locations  in  the  collection  area.  Through  coordination  with  Douglas 
Bickel  at  Sandia  National  Faboratory,  the  laboratory  that  designed  the  Fynx  II  radar,  the 
optimum  radar  cross  section  required  for  best  radar  response  was  determined  to  be 
25dBsm  and  field-portable  trihedral  reflectors  were  designed  and  fabricated,  Figure  13. 
These  field  portable  corner  reflectors  facilitated  the  ability  to  transport  a  total  of  five 
reflectors  and  tripods  in  a  single  backpack  to  the  collection  sites. 
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Figure  13.  Lightweight  portable  corner  reflectors  were  designed  to  be  broken  down  for  use  in 
the  field.  Five  of  these  along  with  their  stands  could  easily  be  put  in  to  a 
backpack  and  either  skied  or  hiked  in  to  the  SAR  collection  areas. 


Two  of  the  eight  identified  sites  were  chosen  for  SAR  collection.  These  sites  will 
be  referred  to  as  “Elysian  Fields”  and  the  “CRREL  snow  study  site.”  Ultimately,  while 
both  sites  were  prepared  for  the  collection,  due  to  flight  time  restrictions,  SAR  data  were 
collected  only  on  Elysian  Fields.  Site  preparation  for  the  Snow  On  collection  took  place 
on  the  first  and  second  of  April  2012.  The  first  step  in  preparing  the  sites  was  to  identify 
the  center  of  the  SAR  collection  areas.  These  coordinates  were  taken  on  the  earlier 
performed  site  survey.  With  these  coordinates  in  mind,  locations  for  the  corner  reflectors 
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were  chosen.  The  initial  plan  for  the  Elysian  Fields  collection  area  was  to  place  the 
comer  reflectors  100  meters  apart  in  the  shape  of  a  square  centered  on  the  target  location 
coordinates.  This  had  to  be  adjusted  during  site  preparation  because  a  100m  square 
encroached  into  the  tree  line  on  the  south  side  and  up  against  the  ski  area  boundaries  on 
the  north  side.  Therefore,  the  reflectors  were  adjusted  to  be  approximately  80  meters 
apart  in  the  north-south  direction  while  the  east-west  direction  remained  at  100  meters 
(Figure  14).  The  snow  was  excavated  to  the  ground  at  each  corner  reflector  location  to 
emplace  the  tripods.  The  GPS  locations  were  recorded  and  a  stake  was  placed  to  mark 
the  location.  The  stake  was  critical  to  allow  for  return  of  the  reflectors  to  the  same 
location  during  the  Snow  Off  SAR  collection.  The  snow  pits  dug  for  the  comer  reflectors 
varied  between  1.5  to  2.3  meters  in  snow  depth  at  this  collection  site.  Three  of  the  four 
tripods  were  positioned  directly  on  the  bare  ground  and  allowed  the  reflector  to  be  at  or 
above  the  snow  line,  making  them  visible  to  the  radar  (Figure  15).  The  fourth  tripod,  in 
the  deepest  pit,  was  positioned  with  its  legs  placed  into  the  side  walls  of  the  pit  at  a  height 
that  allowed  the  reflector  to  be  visible  to  the  aircraft.  All  the  reflector  stands  were  then 
staked  and  tied  down.  The  reflectors  were  attached  to  the  heads  of  the  tripods  on  the  day 
of  the  collection  with  a  clear  view  of  the  sky  towards  the  flight  track.  They  were  pointed 
due  south  or  167°  magnetic  and  with  an  elevation  angle  of  -13°  to  correspond  with  GAA 
SAR  parameters  and  planned  flight-lines.  The  direction  was  detennined  with  a  magnetic 
compass  and  verified  with  a  GPS  compass.  The  elevation  angle  was  accomplished  with  a 
carpenter’s  protractor.  This  proved  to  provide  a  very  well  pointed  reflector  and  the 
aircrew  could  clearly  identify  the  reflectors;  thus  no  additional  adjustments  were 
required. 
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Figure  14.  Corner  reflector  orientation  at  Elysian  Fields  collection  site  during  Snow  On  SAR 
collection  (37°37.7’N,  1 19°02.7’N).  The  image  area  is  312  x  176m  with  the  N-S 
corner  reflectors  80m  apart  and  the  E-W  comer  reflectors  100m  apart.  Magnitude 
image  provided  by  General  Atomics  Aeronautical. 
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Figure  15.  Deployment  of  corner  reflectors  used  as  ground  control  points  at  Elysian  Fields 

collection  site. 


The  Snow  Off  collection  was  similar  to  the  Snow  On  collect.  The  majority  of  the 
site  preparation  was  done  the  day  of  the  collection,  Friday  13  July  2012.  Initially,  the 
trihedral  reflectors  were  set  up  at  the  same  four  locations  as  during  the  Snow  On 
collection.  They  were  mounted  on  tripods  as  before  and  raised  to  the  same  recorded 
elevation  above  the  bare  ground  as  the  initial  collection.  One  tripod,  the  northwest 
corner,  had  to  be  modified.  This  location  is  where  the  tripod  could  not  be  place  on  the 
ground  due  to  snow  depth.  For  this  collection,  an  extension  was  manufactured  for  the 
tripod  legs  out  of  three-sided  channel  aluminum.  This  proved  to  be  stable  and  strong 
enough  to  meet  the  requirements.  Once  the  tripods  were  positioned  and  brought  to  the 
correct  elevation  above  the  ground,  the  reflectors  were  attached  and  pointed  as  before. 
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Five  additional  reflectors  were  put  in  place  during  the  Snow  Off  collect  to  help 
assess  SAR  accuracy  (Figure  17).  They  were  positioned  directly  on  the  ground  and 
pointed  in  accordance  with  the  previous  mentioned  procedures.  Four  of  the  reflectors 
were  placed  approximate  35  meters  out  in  a  diagonal  direction  from  the  four  already 
established  reflectors.  The  fifth  one  was  placed  in  the  center  of  the  square  formed  by  the 
already  established  reflectors.  The  resulting  shape  was  an  “X”  pattern.  It  should  be 
mentioned  that  the  additional  reflectors  were  intended  to  help  improve  the  accuracy  of 
the  interferogram  processing  and  the  subsequent  DEM.  They  were  not  however  directly 
used  as  control  points  when  comparing  the  Snow  On  and  Snow  Off  images  since  they 
were  not  present  during  the  Snow  On  collection. 


Figure  16.  Corner  reflector  deployment  for  Snow  Off  SAR  collection. 
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Figure  17.  Comer  reflector  orientation  at  Elysian  Fields  (3 12  x  176m)  collection  site  during 
Snow  Off  SAR  collection  (37°37.7’N,  1 19°02.7’N).  The  ninth  comer  reflector  is 
outside  the  image  area.  Magnitude  image  provided  by  General  Atomics 

Aeronautical. 
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3.  Site  Survey 

The  process  of  performing  the  calculations  to  detennine  the  baseline  in  the 
Baseline  method  and  the  linear  regression  in  the  BFPR  method  required  as  much 
precision  as  possible  for  the  comer  reflector  locations.  According  to  GPS.Gov  2013,  an 
accuracy  of  three  meters  can  be  expected  when  using  standard  GPS  units.  This  is  under 
optimal  conditions  and  errors  can  easily  be  greater.  To  improve  on  this,  a  survey  grade 
GPS  was  used  to  determine  CR  locations.  The  Trimble  GeoExplorer  6000  (centimeter) 
was  the  specific  unit  used.  Through  post  processing  this  unit  is  advertised  to  have  an 
accuracy  of  1cm  in  the  horizontal  and  1.5cm  vertical  (Trimble  2013).  The  corner 
reflector  locations  were  surveyed  with  this  unit  multiple  times  at  different  times  of  the 
day  for  Snow  Off  conditions.  Post  processing  on  the  surveys  showed  positioning 
consistency  within  a  centimeter  or  two.  Due  to  the  limited  availability  of  this  equipment, 
it  was  only  used  in  the  surveying  of  the  corner  reflector  locations. 

4.  Validation 

The  fourth  aspect  of  the  field  work  was  the  validation.  To  ensure  that  the  SAR 
snow  depth  retrieval  techniques’  accuracy  could  be  tested,  in-situ  snow  depth 
measurements  were  taken  concurrently  with  the  Snow  On  data  collection  on  3  April 
2012.  These  measurements  were  taken  using  a  1cm  graduated  avalanche  probe 
immediately  after  the  General  Atomics  aircrew  finished  collecting  the  datasets  of  the 
Elysian  Fields  area  during  the  Snow  On  collection.  Due  to  the  disturbed  snow  directly 
between  the  corner  reflectors,  the  measurements  were  taken  approximately  every  20 
meters  within  the  box  outlined  by  the  reflectors.  A  total  of  16  measurements  were  taken 
in  a  roughly  80x80m  grid.  Due  to  the  shortening  of  the  box  outlined  by  the  comer 
reflectors,  the  in-situ  snow  depth  grid  extended  outside  of  the  planned  analysis  area 
marked  by  the  corner  reflectors  to  the  south.  In  addition  to  these  measurements,  random 
measurements  were  taken  further  south  of  the  planned  verification  area.  These 
measurements  extended  into  the  tree-line  and  may  provide  some  insight  to  the  impacts  of 
foliage  on  this  method.  All  the  snow  depth  measurement  locations  were  recorded  with  a 
GPS  along  with  calculated  GPS  error.  In  addition  to  these  snow  depth  measurements,  a 
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snow  pit  was  dug  concurrently  with  the  SAR  collection  for  snow  analysis  near  the 
northeast  CR.  It  included  temperature  readings,  crystal  size  and  type,  and  density 
measurements  throughout  the  column. 


Figure  18.  An  avalanche  snow  analysis  pit  was  dug  during  the  SAR  collection.  A  vertical 
profde  of  the  snowpack’s  temperature,  crystal  size,  crystal  type,  and  snow  density 

was  recorded. 


C.  DATA  PROCESSING  (BASELINE  METHOD) 

The  Baseline  InSAR  data  processing  method  is  made  up  of  several  well 
documented  steps.  The  path  used  for  this  method  modifies  the  path  that  is  described  by 
M.  A.  Richards  2007  and  J.  A.  Richards  2009  and  was  carried  out  using  the  MATLAB 
programming  language.  The  following  steps  were  used  in  the  process: 

•  Data  preparation 

•  Complex  image  generation 

•  Image  registration 
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•  Interferogram  generation 

•  Phase  unwrapping 

•  Base  line  detennination 

•  Flat  Earth  removal 

•  Height  to  phase  scaling 

•  DEM  generation 

•  Snow  On/Snow  Off  differencing 

•  Verification 

1.  Data  Preparation 

The  SAR  data,  as  mentioned  earlier,  were  collected  by  General  Atomics.  The 
data  are  recorded  in  proprietary  “.limg”  files  that  can  only  be  read  through  the  General 
Atomics  CLAW  III  imaging  software  (D.  Bickel,  2012,  personal  communication)  .  The 
software  was  not  directly  available  for  this  research  so  Sandia  National  Labs  and  General 
Atomics  pulled  the  relevant  data  out  from  the  Snow  On  and  Snow  Off  collections 
respectively  and  shipped  it  to  NPS  as  MATLAB  files.  The  data  consisted  of  three  pieces 
of  infonnation  for  each  pass.  These  were  the  raw  magnitude  of  the  reflectivity,  the 
recorded  phase,  and  the  header  data  consisting  of  the  specific  flight  and  radar  pointing 
infonnation. 


2.  Complex  Image  Generation 


Once  the  data  were  available  in  a  fonn  that  could  be  used  in  MATLAB,  the 
magnitude  and  phase  arrays  were  combined  to  create  a  complex  image  for  each  pass. 
These  complex  images  are  required  for  the  complex  coherence  calculations  and  are  the 
basis  for  the  generation  of  the  interferogram  (Richards  2009).  The  complex  image  array 
can  be  calculated  as  shown  in  (18),  “j”  represents  the  imaginary  number  square  root  of 

CC  59 


Complex  _  image  =  magnitude  x  e 


( jx  phase)* 


2n 

2“ 


(18) 
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It  should  be  noted  here  that  the  27t/2 16  is  a  conversion  factor  to  account  for  the  fact 
that  the  raw  phase  is  encoded  as  a  16-bit  integer  during  the  SAR  collection  (D.  Bickel, 
2012,  personal  communication). 

3.  Registration 

As  described  in  the  background  section,  the  generation  of  DEMs  from  SAR  data 
requires  two  different  viewing  angles.  For  this  research,  this  was  accomplished  with 
SAR  acquisitions  from  two  separate  individual  passes.  The  arrays  or  images  from  the 
first  pass  will  be  referred  to  as  the  “master”  array/image  and  the  data  from  the  second 
pass  will  be  referred  to  as  the  “slave”  array/image.  Because  the  two  images  are  viewed 
from  different  viewing  angles  they  have  different  perspectives.  They  therefore  must  go 
through  a  registration  process  that  lines  up  the  pixels  from  the  slave  image  with  those  of 
the  master  image  so  that  the  proper  calculations  can  be  made. 

The  first  step  in  the  registration  of  these  arrays  begins  with  identifying  control 
points  in  each  of  the  images.  Because  the  control  points  are  common  between  the 
images,  they  can  be  used  to  help  align  those  images.  As  stated  earlier,  control  points 
were  marked  on  the  ground  using  trihedral  comer  reflectors.  Locating  these  control 
points  in  the  images  was  therefore  a  fairly  straightforward  process.  The  magnitude  array 
was  displayed  as  an  image  showing  the  collection  area.  Once  this  image  was  displayed, 
the  corner  reflectors  were  clearly  visible.  The  values  of  the  pixels  associated  with  the 
center  of  the  reflectors  were  several  orders  of  magnitude  greater  than  those  of  the 
surrounding  pixels.  Being  able  to  clearly  identify  the  location  of  the  reflector  center 
allowed  for  the  identification  of  the  exact  X-Y  location  in  the  array.  These  identified 
locations  were  ultimately  used  in  the  processing  that  warped  the  slave  image  to  the 
master  image,  which  ensures  that  the  pixel  locations  between  the  master  and  slave  images 
line  up  correctly.  This  is  necessary  for  the  interferometric  phase  generation  and  any  other 
direct  comparisons  made  between  images.  The  corner  reflectors  can  be  seen  the  “bright 
stars”  in  the  zoomed  in  image  Figure  19. 
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Figure  19.  Magnified  Snow  On  SAR  magnitude  image  showing  two  comer  reflectors.  In 
addition  to  the  reflectors,  the  walking  paths  can  be  seen  in  the  snow  as  bright 
linear  tracks.  Image  size  45x77m,  Mammoth  Mountain  study  site,  (37°37.7’N, 

1 19°02.7’N). 

The  warping  process  itself  consists  of  four  steps  and  is  achieved  using  several 
specialized  functions  in  the  MATLAB  image  toolbox  (Mathworks,  2013).  Using  the 
identified  control  points  in  both  the  master  and  slave  image  along  with  the  image  toolbox 
“cporr”  function  which  uses  normalized  cross-correlation  to  adjust  each  pair  of  control 
points,  the  slave  control  points  are  tuned  to  the  master  control  points.  This  provides  an 
output  that  is  used  to  generate  a  spatial  transform  that  is  used  to  warp  the  rest  of  the 
image.  The  image  toolbox  “cp2form”  function  uses  the  above  output  to  create  this  spatial 
transform. 
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The  spatial  transform  of  the  slave  image  unfortunately,  as  would  be  expected, 
causes  a  change  in  the  array  size.  This  change  in  array  size  prevents  any  further  array 
calculations.  Therefore  it  is  necessary  that  the  proper  array  size  is  maintained  for 
perfonning  calculations.  This  is  easily  enough  circumvented  using  the  “size”  function  to 
match  the  array  sizes. 

With  the  arrays  the  same  size  and  the  spatial  transformation  generated,  the  last 
step  is  to  perform  the  final  warping  of  the  slave  image.  This  is  accomplished  using  the 
MATLAB  image  toolbox  “imtransfonn”  function,  which  takes  the  previously  calculated 
spatial  transform  and  applies  it  to  the  entire  image. 

4.  Interferogram  Generation 

Once  the  slave  image  is  properly  warped  to  the  master  image,  the  next  step  is  to 
generate  the  interferogram  and  complex  coherence.  The  interferogram  “*”  can  be 
calculated  according  to  (19)  (Richards  2009) 

i(x,y)  =  el(x,y)e2(x,y)  (19) 

The  complex  master  array  and  complex  conjugate  slave  array  are  represented  by 

sfe  9  >  <  •••  • 

“ei”  and  “e2  ”  respectively.  Again,  this  calculation  is  straightforward  in  MATLAB.  An 
example  of  an  interferogram  can  be  seen  if  Figure  20. 
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Snow  Off  21/22  Interferogram 


Image  pixel  position  "x" 


Figure  20.  Example  of  interferogram  generated  during  the  Snow  Off  SAR  collection.  The 
repeating  fringe  pattern  is  derived  from  the  subtraction  of  the  phases  of  the  two 
interferometric  datasets.  Both  datasets  record  values  between  -n  and  n. 
Therefore  the  phase  difference  between  the  fringes  can  never  be  greater  than  27t. 
Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


Displaying  the  interferogram  demonstrates  how  the  terrain  changes  based  on  the 
measured  change  of  phase  from  one  pixel  to  the  next.  Because  phase  difference  is  being 
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measured,  only  values  between  -n  to  n  radians  are  recorded.  The  resulting  “wrapped” 
phase  image  displays  repeating  values  between  -n  to  n  as  was  seen  in  Figure  20. 


In  addition  to  the  interferogram,  the  complex  coherence  is  a  useful  parameter  to 
calculate  as  it  provides  infonnation  about  the  quality  of  the  radar  returns  and  their 
suitability  for  generating  DEMs.  The  complex  coherence  represented  by  “y”  can  be 
calculated  as  seen  in  (20)  (Richards  2009). 


An  example  of  complex  coherence  can  be  seen  in  Figure  2 1 . 


45 


Coherence  21/22 


1 


Figure  2 1 .  Example  of  coherence  image  generated  from  the  Snow  Off  SAR  collection. 

Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


5.  Phase  Unwrapping 

The  “wrapped”  interferometric  phase  image  needs  to  be  “unwrapped”  to 

detennine  the  overall  phase  difference  between  the  two  SAR  images  and  to  generate  an 

image  displaying  constant  change  in  phase.  The  phase  unwrapping  process  consists  of 

changing  the  2rc  cyclical  variations  in  phase  to  a  constant  change  in  phase  that  will 

eventually  be  used  to  detennine  elevation  change  across  the  scene  (Richards  2009). 
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Unwrapping  the  interferogram  generated  from  the  SAR  image  pair  produces  a 
single  array  or  image  showing  the  change  in  phase  attributable  to  terrain.  It  is  important 
to  note  that  because  this  array  is  measured  in  phase,  that  the  units  being  used  to  describe 
the  change  in  the  terrain  are  radians.  The  branch  cut  method  developed  by  R.M. 
Goldstein  et  al.  1988  was  used  to  unwrap  the  Mammoth  Mountain  interferograms.  The 
MATLAB  unwrapping  algorithms  implementing  this  approach  were  written  by  Bruce 
Spottiswoode  and  latter  modified  by  Carey  Smith  (Mathworks  2013)  to  increase 
processing  efficiency.  An  example  of  the  unwrapped  phase  from  Figure  20  can  be  seen 
in  Figure  22. 
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Figure  22.  Unwrapped  total  phase  generated  from  interferometric  image  collected  during 
Snow  On  SAR  collection.  The  repeating  2jt  pattern  has  been  removed  and  the  full 
range  of  phase  difference  can  be  seen.  Mammoth  Mountain  study  site, 
(37°37.7’N,  119°02.7’N). 


The  above  interferometric  image  is  not  yet  ready  to  be  converted  to  elevation.  It 
consists  of  phase  that  is  due  to  changes  in  the  x,  y,  and  z  directions.  The  changes  in  the 
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phase  due  to  changes  in  the  x  and  y  direction  is  often  referred  to  as  the  “Flat  Earth 
Phase”.  The  total  phase,  which  is  represented  by  the  unwrapped  phase  in  Figure  20,  is  a 
combination  of  the  phase  change  from  the  elevation  change  that  we  are  trying  to  derive 
and  the  flat  earth  phase,  which  can  be  shown  in  (21),  where  “A0totai”  is  the  total  phase 
seen  in  Figure  22.  “A 0Z”  is  the  phase  change  associated  with  the  change  in  the  elevation, 
and  “A <pfiat_earth”  is  the  flat  earth  phase  change  (Richards  2009).  For  simplicity, 
atmospheric,  noise,  and  baseline  uncertainty  effects  on  phase  are  not  considered  at  this 
time. 

^ 'total  =  +  'jjat  _earth  ^ 

Therefore,  to  get  at  the  elevation  at  each  pixel,  the  flat  earth  phase  must  be 
removed.  To  do  this,  the  flat  earth  phase,  which  requires  precise  knowledge  of  the  base 
line,  must  first  be  calculated,  and  is  the  subject  of  the  next  sub-section. 


6.  Baseline  Determination 


Considering  the  interferometric  phase  equation  (equation  10,  repeated  below),  it 
can  be  seen  that  the  interferometric  phase  is  highly  sensitive  to  the  baseline  length  “B”. 
The  baseline  is  used  in  the  calculation  of  the  flat  earth  phase  that  will  be  removed  from 
the  total  phase  and  the  height  to  phase  ratio  to  get  to  the  finale  DEM.  According  to 
Richards  2007,  the  baseline  length  needs  to  be  known  to  an  accuracy  of  a  factor  of  10’4  to 


10'6  of  the  absolute  slant  range  to  the  target. 


A</> 


AnB  sin  0 

I 


(10) 


The  level  of  accuracy  required  for  the  baseline  can  be  a  bit  problematic.  If  the 
SAR  images  are  being  collected  from  ranges  on  the  order  of  10,000  meters,  the  baseline 
accuracy  needs  to  be  sub-meter  at  a  minimum  with  cm  accuracy  preferred.  Aircraft 
positioning  is  generally  handled  by  GPS.  Recall  that  standard  GPS  accuracy  does  not 
approach  the  level  of  precision  that  is  needed  for  the  baseline  determination.  Consider 
the  fact  that  the  aircraft’s  position  is  being  determined  for  both  passes  at  an  accuracy  of 
around  three  meters,  it  is  possible  that  the  GPS  error  can  be  on  the  order  of  six  meters 
under  optimal  conditions.  Because  of  this,  there  needs  to  be  a  better  way  to  estimate  the 
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baseline  between  the  flights.  It  may  be  possible  to  use  the  SAR  data  themselves  to 
determine  the  baseline  if  there  is  infonnation  about  the  three  dimensional  change 
available  for  multiple  locations  in  the  target  area. 


To  begin  solving  for  the  baseline,  go  back  to  the  interferometric  phase  equation 
which  can  be  seen  in  expanded  form  in  (22). 


.  ,  An  .  .  An  A  AnB  sin  0 

A</>  =  ~  fo  “  ri  )  =  ~YAr  =  J 


(22) 


Recall  that  the  goal  in  the  derivation  of  the  interferometric  phase  angle  is  to 
determine  the  difference  in  the  length  between  the  two  ranges  as  it  is  measured  in  radians. 
The  earlier  simplified  derivation  must  be  expanded  to  consider  the  fact  that  in  airborne 
SAR  collections,  there  is  no  way  to  insure  that  the  two  collections  are  taken  in  a  flat  plane 
relative  to  the  earth  as  was  shown  in  the  earlier  example.  The  two  flight  paths  orientation 
is  more  likely  to  be  something  along  the  lines  seen  in  Figure  23. 
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Figure  23.  Geometry  for  Multi-pass  Single  Look  Complex  interferometric  SAR  collection. 

“ri”  and  “r2”  are  the  respective  ranges  from  antennas  1  and  2.  “B”  represents  the 
baseline  between  the  two  antenna  locations.  “BA  represents  the  orthogonal 
baseline  between  the  two  radar  beam  paths.  SRP  is  the  scene  reference  point.  “9” 
represent  the  angle  between  the  baseline  and  the  orthogonal  baseline.  “5a” 
represents  the  angle  between  the  two  range  vectors.  (After  Richards  2009;  with 
input  from  Douglas  Bickel  personal  communication,  2013). 


Due  to  the  SAR  geometry;  there  are  multiple  components  of  the  phase  that  need 
to  be  considered.  Keep  in  mind  that  phase  in  the  above  figure  is  merely  a  rendition  of  the 
difference  in  range  between  the  two  different  collections.  The  radar  range  from  one  pixel 
to  the  next  is  determined  by  three  different  aspects.  It  is  determined  by  the  difference  in 
the  change  in  the  range  in  the  x,  y,  and  z  directions.  The  previous  phase  image  consists  of 
phase  that  is  made  up  of  changes  in  these  three  directions  and  the  change  in  the  overall 
range  in  3-dimensional  space  will  be  different  between  the  two  aperture  centers  from  one 
pixel  to  the  next.  That  is  good  because  it  allows  for  the  generation  of  a  3-dimensional 
image.  That  also  means  though,  that  if  extraction  of  elevation  information  is  the  goal,  the 
impact  that  the  change  in  the  x  and  y  directions  will  have  on  the  interferometric  phase, 
which  is  the  aforementioned  flat  earth  phase  will  need  to  be  considered.  Looking  at 
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Figure  8  again,  notice  the  changes  in  either  the  y  and  z  direction  on  the  ground  affect  the 
ranges  to  the  two  aperture  centers  which  ultimately  affect  the  interferometric  phase. 


Recall  from  equations  7  through  9  that 

i\  -  r2  cos  SO  +  B  sin  9 

Using  the  small  angle  approximation 

t\  =r2+  B  sin  0 


Therefore 


(23) 

(24) 


A r  =  r]  —r2  —  B  sin  0  (25) 

Also,  recall  From  Figure  23  that  B  sin(9)  is  merely  the  length  of  the  portion  of 
“rf’  that  goes  from  the  intersection  of  the  perpendicular  baseline  to  “Aperture  Center  1”. 
Notice  that  this  length  can  also  be  characterized  as  the  projection  of  vector  B  onto  vector 
r1.  This  can  be  calculated  by  taking  the  dot  product  of  vector  B  onto  the  unit  vector  of?! 
or  Uri. 


B  sinO  —  B  ■  UVl 

Where 

(26) 

B  —<  Bx,By>Bz,> 

And 

(27) 

The  magnitude  of  rx  can  be  given  by 

(28) 

r  =  Vu2+U2  +rzz 

Resulting  in 

- *  vZ  Tv  rZ 

(29) 

U  -<  1£  JL  i£  > 

1  r  r  r 

Which  means  that 

(30) 

jj  * _  Bxrx+Byry+Bzrz 

ri  r 

(31) 

Making  a  substitution,  the  interferometric  phase  angle  can  now  be  defined  as: 


A0 


471  Bxrx+Byr-y+Bzrz 

7  r  ~ " 


(32) 
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The  equation  for  the  interferometric  phase  angle  at  a  given  location  (x,  y,  z)  in 
terms  of  that  location  and  the  components  of  the  baseline  has  now  been  defined. 
Unfortunately  there  is  still  a  lot  of  infonnation  needed.  A(j),  X,  and  r  are  known  but  Bx, 
By,  and  Bz  are  unknowns. 


z 


Figure  24.  Comparing  the  physical  change  in  location  and  phase  of  multiple  points  to  a 
reference  point  “RP”  enables  the  construction  of  a  system  of  equations  to  solve 
for  the  baseline.  The  given  ranges  to  the  reference  point  “ro”  and  the  other  points 
“ra”,  is  also  required  for  the  solving  the  system  of  equations.  (After  Douglas 
Bickel,  personal  communication,  2013. 


If  a  second  point  “a”  is  chosen  at  location  (xa,  ya,  za)  as  seen  in  Figure  24,  “a”  can 
be  written  as 

a{x,y,z)  =  (x0  +  Axa,  y0  +  ya,  z0  +  Aza)  (33) 

Therefore,  when  Bx,  By,  Bz,  represent  the  components  of  the  baseline;  x0,  yo,  z0, 
are  the  coordinates  of  the  reference  point;  Axa,  Aya,  Aza,  are  the  x,  y,  z  distance  to  point 
“a”;  and  “r”  in  the  range  to  target;  Atj)  at  point  “a”  becomes 
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(34) 


A^fl  = 


4;r  Bx (x0  +  Ax a)  +  By (y0  +  Aya )  +  Bz(z0  +  Aza ) 
A  r 


Separating  the  components  of  (34)  results  in 


A  /  A  Bxx0  +By0+Bzz0  BxAxa  +BvAya  +B:Aza 

—  = - = - + - ^ - 

4  nr  r  ft  (35) 

Focusing  on  the  first  component  of  the  right  hand  side  of  (35),  notice  that  the 

variables  are  unknown,  but  it  is  known  that  this  is  a  constant  and  that  it  is  equal  to  the 

value  of  at  the  reference  point  (RP)  (xo,  yo,  zo)  since  A(x,y,z)  is  0.  The  second 


component  of  the  right  hand  side  of  (35)  is  responsible  for  all  of  the  change  in  the  phase 
relative  to  the  RP.  Therefore  if  Acj)  is  known  at  the  RP  and  three  different  points  along 
with  their  relative  A(x,y,z),  there  would  be  three  equations  and  three  unknown  variables 
that  allow  for  a  system  of  equations  to  solve  for  those  variables.  This  approach  to  solving 
for  the  baseline  was  developed  in  cooperation  with  Douglas  Bickel  of  Sandia  National 
labs.  Figure  25  depicts  the  four  locations  being  used  in  the  system  of  equations. 
MATLAB  algorithms  designed  to  solve  multiple  simultaneous  equations  were  used  to 
make  these  calculations. 
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Figure  25.  This  approximately  200x300m  scene  depicts  a  reference  point  ”RP”  defined  by 
the  SE  comer  reflector  and  the  relative  positions  of  the  other  three  known 
positional  points  defined  by  the  other  three  emplaced  comer  reflectors  at  positions 
“a”,  “b”,  and  “c”.  Mammoth  Mountain  study  site  (312  x  176m),  (37°37.7’N, 

1 19°02.7’N). 
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Before  the  system  of  equations  can  be  solved,  the  data  used  to  solve  these 
equations  needs  to  be  addressed.  The  interferometric  phase  angle  is  straightforward.  A(j) 
has  already  been  calculated  as  discussed  in  previous  sections.  The  real  world  A(x,y,z) 
also  needs  to  be  detennined.  As  discussed  earlier,  standard  GPS  does  not  provide  the 
level  of  accuracy  needed  to  make  these  calculations.  Therefore  the  survey  grade  GPS  was 
used  to  determine  the  A(x,y,z). 

Knowing  the  real  world  A(x,y,z),  is  not  enough  though.  The  real  world  “y” 
coordinate  orientation  is  true  north/south  and  the  real  world  “x”  coordinate  orientation 
true  east/west.  Despite  best  efforts  and  potential  other  collection  scenarios;  the  SAR 
image  coordinate  system  does  not  match  the  real  world  coordinate  system.  This  can  be 
seen  in  Figure  26.  Because  all  calculations  are  being  performed  in  a  SAR  image  centric 
coordinate  system,  the  real  world  A(x,y)  needs  to  have  its  coordinate  system  rotated  to 
match  that  of  the  SAR  acquisition.  When  this  is  done,  the  Ax  and  Ay  values  will  change 
and  will  therefore  need  to  be  recalculated  based  on  the  coordinate  system  rotation. 
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Radar  image  area  surface  projection 


N(true)  \ 


Figure  26.  In  order  to  correctly  determine  the  baseline,  an  accurate  representation  of  the  x,  y, 
and  z  changes  must  be  known.  Those  changes  were  computed  in  a  reference 
frame  where  the  x  axis  was  east/west  and  the  y  axis  was  north/south.  The  radar 
reference  frame  does  not  line  up,  thereby  requiring  an  adjustment  to  the  x  and  y 
components  in  the  baseline  calculation  equations. 


Reorientation  from  the  real  world  coordinate  system  to  the  SAR  coordinate 
system  is  a  matter  of  geometry/trigonometry  and  can  be  seen  in  Figure  27. 
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p[  =  pl  cos(ci')  +  p2  sin(or)  p\  =  -px  sin(ar)  +  p~,  cos(or) 

Figure  27.  Demonstration  of  the  correction  required  for  a  given  rotation  of  a°  rotation  in 
reference  frames  (From  Mathematics-Online  Lexicon,  2013). 

To  determine  the  proper  rotation  to  calculate  the  new  A(x,y),  “a”  must  first  be 
determined  in  the  above  equations.  That  is  where  the  third  piece  of  SAR  data  provided 
by  Sandia  and  General  Atomics  comes  in.  Recall  that  third  data  set  for  each  SAR 
acquisition  was  the  “Header  Data”  that  provides  both  flight  information  and  radar 
pointing  information.  An  example  of  the  header  data  can  be  seen  in  Figure  28. 
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Header  data:  Scan  21  0.1m 


Field  ^ 

Value 

Min 

Max 

dd  version 

513 

513 

513 

dd  submode 

7 

7 

7 

di  bytes_per_pixel 

4 

4 

4 

n  navtime 

4.9672e+05 

4.9672... 

4.9672... 

jid  r_b 

1.0293e*04 

1.0293... 

1.0293... 

iB  x_dist 

1.1531e+03 

1.1531... 

1.1531... 

dd  y_dist 

-9.9595e+03 

-9.959... 

-9359... 

n  z_dist 

-2.3295e+03 

-2.329... 

-2.329... 

Iv.x 

90.1654 

90.1654  90.1654 

S  V_y 

0.7233 

0.7233 

0.7233 

dd  v_z 

0.0956 

0.0956 

0.0956 

dd  vlos 

10^226 

103226  103226 

dd  phi 

13.0802 

13.0802  13.0802 

_U  cos_phi 

0.9741 

0.9741 

03741 

dd  theta 

82.9342 

823342  82.9342 

dd  lat 

37.5398 

37.5398  37.5398 

±i  long 

-119.0328 

-119.0... 

-119.0... 

S  hgt 

5.4263e+03 

5.4263... 

5.4263... 

dd  heading 

-89.1830 

-89.18... 

-89.18... 

dd  ground_speed 

90.1904 

90.1904  90.1904 

dd  v_north 

1.2860 

1.2860 

1.2860 

dd  v_east 

-90.1812 

-90.18... 

-90.18... 

dd  v_down 

0.1013 

0.1013 

01013 

dd  week 

1696 

1696 

1696 

dd  azpixels 

2048 

2048 

2048 

dd  rngpixels 

3642 

3642 

3642 

dd  azpixspace 

0.0858 

0.0858 

0.0858 

dd  rngpixspace 

0.0858 

0.0858 

0.0858 

dd  rhoa 

0.1016 

0.1016 

0.1016 

dd  rhor 

0.1016 

0.1016 

0.1016 

dd  calparm 

4.8928e-04 

4.8928... 

4.8928... 

dd  sigO_eff_db 

-35.7972 

-35.79... 

-35.79... 

dd  scla 

263.7491 

263.74... 

263.74... 

dd  range_nom 

1.0302e«-04 

1.0302... 

1.0302... 

i  j  phi_nom 

13.1147 

13.1147  13.1147 

dd  theta_nom 

833976 

83.3976  83.3976 

dd  ap_count 

ill  srp 

0 

<1x1  struct> 

0 

0 

Figure  28.  Header  data  provided  with  each  SAR  acquisition  provides  specifics  on  the  aircraft 
flight  path  orientation  and  the  orientation  of  the  radar  with  respect  to  the  flight 
information.  In  addition  to  the  flight  data,  the  header  provides  the  calculated 
pixel  size  for  the  given  image  or  the  ground  resolution. 
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There  are  two  key  pieces  of  information  with  respect  to  determining  “a”.  The 
first  number  of  interest  is  the  “sola”  which  represents  the  flight  path  of  the  aircraft 
relative  to  true  north.  In  the  above  header  data,  “sola”  can  be  seen  to  be  263.7491°.  The 
“sola”  can  be  visualized  in  Figure  29. 


N 


A 


Figure  29.  SCLA  represents  the  flight  path  relative  to  true  north. 

The  second  piece  of  information  that  is  required  is  “theta.”  “Theta”  represents  the 
squint  angle  or  the  pointing  direction  of  the  radar  relative  to  the  flight  path.  In  the  above 
header  data  “theta”  can  be  seen  to  be  82.9342°.  This  particular  theta  should  not  be 
confused  with  the  use  of  9  in  other  portions  of  this  method.  Due  to  use  of  it  in  the  header 
data  by  General  Atomics  it  will  simple  be  referred  to  as  “theta”  rather  than  the  symbol  in 
the  text  and  will  not  be  renamed.  Theta  can  be  visualized  in  Figure  30. 
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0=  82.9342“ 


Flight  Path  (SCLA=263.7491° ) 


< - - 

Figure  30.  Theta  in  the  header  data  represents  the  squint  angle  or  the  angle  that  the  radar  is 
looking  relative  to  the  flight  path.  In  this  case  theta  is  82.9342°. 

The  next  step  is  to  combine  the  “sola”  and  “theta”  data  to  determine  how  the  SAR 
pointing  direction  deviates  from  true  north.  This  deviation  will  be  the  “a”  used  to 
determine  the  change  in  A(x,y).  The  final  geometry  can  be  seen  if  Figure  3 1  and  the 
equation  used  can  be  seen  in  (36)  resulting  in  a  calculated  a  of  13.3167  in  this  example. 
If  the  flight  path  is  in  the  opposite  direction  and  the  radar  is  pointed  to  the  left  of  the 
direction  of  flight,  “theta”  will  be  a  negative  number  and  (36)  will  have  to  be  adjusted 
accordingly  for  the  differing  geometry. 


oc=  360  —  (scla  +  theta ) 


(36) 
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Figure  3 1 .  This  figure  shows  all  the  relevant  angles  required  to  detennine  the  rotation  angle 
“a”  used  in  calculating  the  change  in  x  and  y  for  the  rotation  of  the  reference 
frame.  “0”  represents  the  squint  angle.  “0(tme)”  represents  the  squint  angle  adjust 
for  true  north.  SCLA  represents  the  flight  heading. 


Now  that  the  correct  “a”  has  been  determined,  A(x,y)  can  be  adjusted  so  that  it 
can  be  used  to  calculate  the  components  of  the  baseline  in  the  system  of  equations  that 
was  discussed  above. 

7.  Flat  Earth  Removal 

Once  the  components  of  the  baseline  have  been  determined,  all  the  information  is 
available  to  determine  the  flat  earth  interferometric  phase.  Recall  that  the  total  phase 
generated  in  the  interferogram  is  a  combination  of  phase  due  to  the  change  in  topography 
and  the  change  in  the  phase  due  to  the  change  in  range  from  the  two  acquisitions  in  the  x 
and  y  directions.  The  flat  earth  phase  needs  to  be  removed  before  the  terrain  phase  can 
be  addressed. 

To  detennine  the  flat  earth  phase,  refer  back  to  the  interferometric  phase  equation. 
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A<f)  = 


4ttB  sin  9 


&  (10) 
The  flat  earth  phase  is  simply  a  representation  of  what  the  interferometric  phase 
would  be  if  the  ground  was  a  flat  plane.  It  can  be  seen  from  (10)  that  this  phase  depends 
on  the  baseline  B  which  is  determined  by  the  flight  geometry  in  3-dimensional  space. 
The  flat  earth  phase  itself  is  made  up  of  two  aspects,  to  include  the  range  direction  and 
the  azimuth  direction.  The  range  aspect  of  the  flat  earth  is  detennined  by  the  height  of 
the  platform  above  the  flat  plane  and  the  ground  distance  in  the  y  direction.  The  azimuth 
component  is  composed  of  strictly  the  x  direction,  the  along  the  flight  path  direction. 


Starting  with  the  range  portion  of  the  flat  earth,  equation  10  shows  that  the 
interferometric  phase  is  dependent  on  three  variables.  The  wavelength  X  is  a  constant  and 
is  known.  The  baseline  “B”  was  calculated  in  the  previous  section  and  is  also  constant 
throughout  the  scene.  The  third  variable  is  the  incidence  angle  9.  Figure  32  shows  that 
when  moving  along  the  established  flat  plane  in  the  y  direction,  that  the  height  remains 
constant  while  9  changes  with  the  change  in  y.  Therefore,  the  knowledge  of  the  change 
in  9  is  required  to  determine  the  change  in  the  interferometric  phase  angle  throughout  the 
scene. 
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Figure  32.  Geometry  shows  that  0  changes  as  the  distance  from  the  antenna  “ro”  changes. 

This  change  in  0  is  required  to  determine  the  interferometric  phase  angle  and  is 
calculated  using  the  known  altitude  “h”  above  the  scene  reference  point  (SRP) 

and  the  ground  distance  “y”. 

The  first  step  in  this  process  is  to  establish  an  array  the  same  size  as  that  of  the 
interferometric  images  being  analyzed.  In  order  to  calculate  0;  both  the  height  and 
ground  range  “y”  must  be  known.  It  should  be  mentioned  at  this  point  that  the  SAR 
system  is  very  good  at  determining  range  to  the  reference  point.  This  allows  accurate 
determination  of  the  platfonn  height  and  the  ground  distance  to  the  center  pixel  in  the 
target  image.  The  center  pixel  in  a  SAR  image  represents  that  target  and  is  referred  to  as 
the  scene  reference  point  (SRP).  All  other  pixels  in  the  image  are  at  known  ranges  but 
the  y  and  z  coordinates  are  not  known  because  both  y  and  z  changes  in  the  terrain 
contribute  to  the  range  at  each  individual  pixel.  The  radar  system  calculates  the  size  of 
the  pixels  in  the  image,  or  the  resolution,  and  provides  it  in  the  header  data.  With  this  in 
mind,  an  imaginary  flat  plane  is  established  through  the  SRP  ensuring  that  the  height  “h” 
is  known.  In  addition  to  the  height,  the  ground  distance  is  known  to  the  reference  point. 
Using  the  ground  distance  to  the  SRP  and  the  processed  ground  resolution,  the  ground 
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distance  to  every  pixel  can  be  determined  in  an  imaginary  plane.  This  allows  calculation 
of  9  for  every  point  on  the  flat  earth  plane.  With  the  calculation  of  9  for  all  of  the 
positions  in  the  array,  the  calculation  of  the  interferometric  phase  with  respect  to  the  y-z 
direction  is  not  difficult.  Simply  apply  the  interferometric  phase  equation  with  X,  the  Byz 
baseline  component,  and  the  newly  calculated  9.  An  example  of  this  phase  can  be  seen  in 
Figure  33  fully  unwrapped  and  Figure  34  wrapped. 
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Figure  33.  This  notional  figure  demonstrates  the  unwrapped  change  in  phase  in  the  “y” 
direction  due  to  the  change  in  range  (range  aspect  of  the  flat  earth  phase). 
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Snow- Off 

Wrapped  y-z  flat  plane  phase 
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Figure  34.  This  notional  figure  demonstrates  the  wrapped  change  in  phase  in  increments  of 
2n  radians  in  the  “y”  direction  due  to  the  change  in  range  (range  aspect  of  the  flat 

earth  phase).  . 

The  second  step  in  determining  the  flat  earth  phase  is  to  determine  the  azimuthal 
phase  change.  The  azimuthal  flat  earth  change  uA<pbxjiat”  can  be  given  by  (37)  where 
“Bx”  is  the  x  component  of  the  baseline  and  “X”  is  the  azimuthal  distance  from  the  SRP. 

_  4  ttBx  X 

(37) 

The  resulting  unwrapped  image  can  be  seen  in  Figure  35  and  the  wrapped  image 
can  be  seen  in  Figure  36. 
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Snow  Off 
Azimuth  Phase 


Figure  35. 
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This  notional  figure  demonstrates  the  unwrapped  change  in  phase  in  the  “x”  or 
azimuth  aspect  of  the  flat  earth  phase. 
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Figure  36.  This  figure  demonstrates  the  unwrapped  change  in  phase  in  increments  of  2n 
radians  in  the  “x”  or  azimuth  aspect  of  the  flat  earth  phase. 

Once  the  flat  earth  phase  in  the  range  (y-z  direction)  and  the  azimuth  direction  are 
determined,  they  can  be  added  together  to  get  the  total  phase  correction  or  the  flat  earth 
correction.  This  can  be  seen  in  Figure37. 
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Image  pixel  position  "x"  Image  pixel  position  "x"  |magc  pjxe|  posltion  "x" 


Figure  37.  The  flat  earth  phase  is  sum  of  the  range  and  azimuth  phases.  Subtracting  the  flat 
earth  phase  from  the  total  unwrapped  phase  results  in  the  terrain  phase. 

The  final  step  is  to  remove  the  total  flat  earth  phase  from  the  total  interferometric 
phase.  Figure  38  (right)  represents  the  interferometric  phase  that  is  due  to  elevation 
change  or  the  topography  alone.  This  is  also  known  as  the  flattened  phase. 


Figure  38.  Subtracting  the  flat  earth  phase  from  the  total  unwrapped  interferometric  phase 
results  in  a  phase  array  or  flattened  phase  that  represents  terrain  phase  only  (right 

image). 
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8.  Height  to  Phase  Scaling 

The  interferometric  phase  due  to  the  terrain  is  the  only  phase  left  once  the  total 
phase  has  been  flattened.  The  next  step  is  to  convert  this  flattened  phase  into  elevation  or 
a  DEM  (Figure  39)  (Rosen  2009). 
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Figure  39.  Geometry  involved  in  determining  phase  to  elevation  scaling  where  “ro” 

represents  the  radar  range  to  an  imaginary  plane,  “T”  represents  a  point  on  the 
terrain  that  is  at  the  same  range  as  ro  but  at  a  height  of  “z”  above  the  imaginary 
plane.  “9o”  represent  the  incidence  angle  to  a  given  point  “0”  on  the  imaginary 
plane  that  is  at  the  same  range  as  that  of  point  T.  “59”  is  the  angle  between  the 
vectors  to  point  O  and  T.  “a”  is  the  angle  between  the  baseline  in  the  y-z  plane 
and  the  horizontal  (After  Rosen  2009;  with  input  from  Douglas  Bickel,  personal 


communication,  2013). 


To  start  the  phase  to  elevation  conversion  process  (Rosen  2009),  Figure  39, 
establishes  the  equations  of  the  imaginary  plane’s  interferometric  phase,  “Atj)o”,  and  the 
total  interferometric  phase  equation,  “A(j)T”,  given  by  (38)  and  (39) 


70 


(38) 


(39) 


Air 

A  (j)T  =  — B  cos(90  -  60  -  86  +  a) 

A 

Expanding  (39)  and  using  the  small  angle  approximation  results  in, 

4  7T 

A(f>T  -  —  B sin(#0  +56 -a) 

A 


(40) 


Two  things  should  be  noted  from  the  above  equations.  The  first  is  that  a  can  be 
detennined  from  the  previously  calculated  baseline  components.  The  second  is  that 
“00+50”  on  the  right  hand  side  of  (40),  is  equal  to  0. 

Expanding  (40)  becomes, 


Aw  Air 

A(f>T  =  — .Ssin^o  -a)  +  — Bcos(60  -a)56 

A  A 


(41) 


Notice  that  the  first  term  in  (40),  is  the  same  as  “A(|>o”.  Therefore,  (41)  can  be 
rewritten  as. 


Rearranging  (42)  results  in, 


4  Tl 

A<j)T  =  A(/>0  + — Bcos(60  -a)56 

A 


4k 


(42) 


A(f>T  -A^0  = — Bcos(60  -a)56  =  A<f>FE 

A 


(43) 


Equation  (43)  represents  the  phase  that  has  had  the  non-terrain  phase  removed  or 
the  same  phase  that  has  been  previously  shown  in  Figure  38. 

To  summarize,  “A(|>t”  is  known;  it  is  the  total  interferometric  phase  calculated  in 
Figure  22.  Likewise  “A(j)0”  is  the  interferometric  phase  calculated  from  the  imaginary 
plane  seen  is  Figure  37.  0o  is  known  from  the  recorded  header  data.  The  baseline  “B”  has 
been  calculated  and  “a”  can  be  detennined  from  the  components  of  “B”.  The  wavelength 
“k”  is  also  known.  The  only  variable  not  known  in  (43)  is  50.  50  can  therefore  be 
calculated  as  seen  in  (44)  and  (45). 

A  A(f>T-A(/)0 

(44) 


SO  = 


4 nB  cos (O0  -  a) 


or 
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(45) 


80  =  — - ^ — 

4 kB  cos (0o  -  a) 

The  final  step  in  deriving  the  phase  to  height  ratio  can  be  seen  in  Figure  40. 


z 


Figure  40.  “50”  or  the  angle  between  the  vectors  that  point  toward  points  “O”  and  “T”  is 

required  to  determine  the  phase  to  height  relationship.  An  equation  can  be 
derived  by  introducing  “n”  for  the  distance  between  “O”  and  “T”  and  recognizing 
that  “9o”,  incidence  angle,  is  the  same  as  the  angle  between  “n”  and  the  imaginary 
plane.  (After  Rosen  2009;  with  input  from  Douglas  Bickel,  personal 

communication,  2013) 


It  can  clearly  be  seen  that, 


sin  89  =  — 

Applying  the  small  angle  approximation,  (46)  becomes, 

89  =  - 

ri 


which  can  be  rearranged  to  be, 

n  =  rx59 


(46) 


(47) 
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(48) 


It  can  also  clearly  be  seen  that, 


sin0o  = 


z 

n 


which  can  be  rearranged  to  be, 


(49) 


z  =  n  sin  60 


(50) 


Combining  (49)  and  (50)  results  in  the  equation  that  provides  the  elevation, 

z  =  r\56  sin  60 

(51)  is  the  final  calculation  needed  to  generate  the  DEM.  Applying  these  calculations  to 
the  previously  interferometric  phase  fields  results  in  Figure  41.  The  mathematical  path  to 
determine  the  phase  to  height  ratio  was  developed  in  cooperation  with  and  primarily  by 
Douglas  Bickel  from  Sandia  National  Labs. 
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Snow  Off  Terrain 
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Figure  41.  Snow  Off  DEM.  Mammoth  Mountain  study  site  (37°37.7’N,  119°02.7’N). 

9.  Snow  On  Snow  Off  Differencing 

The  last  step  in  deriving  the  snow  depth  is  to  subtract  a  DEM  generated 
from  a  Snow  Off  scene  from  a  DEM  generated  from  a  Snow  On  scene.  The 
resulting  image  is  the  result  of  the  change  in  elevation  strictly  due  to  the  addition 
of  snow. 

10.  Baseline  Method  Summary 

Figure  42  is  a  flow  chart  that  summarizes  the  steps  and  processing  order  of 
the  Baseline  method. 
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Baseline  Method 


Figure  42.  Baseline  method  flow  chart. 
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D.  DATA  PROCESSING  (BEST  FIT  PLANE  REMOVAL  METHOD) 

The  “Best  Fit  Plane  Removal  (BFPR)  Method”  takes  a  different  approach  to 
calculating  the  DEM  than  the  Baseline  method.  The  BFPR  approach  does,  however, 
share  some  of  the  same  steps  that  were  performed  in  the  Baseline  method.  The  Baseline 
method  takes  a  direct  route  to  calculating  and  removing  the  flat  earth  phase  “A(|)fiat  earth”, 
and  the  phase  to  elevation  ratio.  These  calculations  require  knowledge  of  the  baseline. 
The  BFPR  method  bypasses  the  baseline  knowledge  requirement,  thus  removing  some 
steps  and  adding  several  others.  The  BFPR  processing  steps,  accomplished  in  MATLAB 
are  as  follows  (shared  steps  italicized): 

•  Data  preparation 

•  Complex  image  generation 

•  Image  registration 

•  Interferogram  generation 

•  Phase  unwrapping 

•  Best  Fit  Plane  generation 

•  Best  Fit  Plane  Removal 

•  Terrain  slope  generation  from  low  resolution  DEM 

•  Control  point  elevation  difference  from  average  slope 

•  Phase  to  elevation  linear  regression 

•  Application  of  linear  regression  equation  to  BFPR  image 

•  Add  linearized  field  to  average  slope 

•  Snow  On/Snow  Off  differencing 

•  Verification 

BFPR  is  best  described  from  a  perturbation  or  decomposition  perspective. 
Variables  with  an  over-bar  represent  the  mean  of  that  variable  while  the  “prime”  symbol 
or  accent  mark  represents  the  deviation  of  the  value  from  a  particular  mean.  In  the 
following  equations  “A^tar  represents  the  total  phase  or  the  unwrapped  interferometric 
phase,  "Ac|)fiat  earth"  represents  the  flat  earth  phase,  "A(|)z”  represents  the  phase  associated 
with  the  terrain. 
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Recall  from  equation  (21)  that  unwrapped  interferogram  is  made  up  of  both  flat 
earth  and  terrain  phase  as  is  reiterated  in  (52)  (Richards  2009). 

Afitotal  At  flat  earth  Ar z  (52) 

From  a  perturbation  perspective  the  flat  earth  phase  is 

Afi flat  _  earth  Afi flat  _  earth  +  ^ flat  _  earth  (53) 

The  flat  earth  phase,  however,  is  a  plane  and  has  no  perturbation.  Therefore  it 
reduces  down  to 

Afi flat  _  earth  A^ flat  _  earth  (54) 

The  terrain  phase  from  a  perturbation  perspective  is 

Af.  =  tfl).  +  A  <j)[  (55) 

Unlike  the  flat  earth  phase,  there  are  variations  throughout  the  image.  “Ac/>z” 
represents  the  average  slope  of  the  terrain  and  “A 0Z'”  is  the  variation  or  perturbation 
from  that  average  slope. 

By  replacing  (54)  and  (55)  into  (53),  the  total  phase  can  now  be  given  as 

A^ total  ~  A^ flat _  earth  +  A  (4  +  A  (4 

or  (56) 

Afi total  —  A ^ flat  _  earth  +  Af4+A<4 

Taking  the  best  fit  plane  of  the  total  phase  is  the  same  as  finding  the  average  slope 
of  the  phase  image  and  is  now  given  by 

Afi total  ~  A^ flat  _  earth  +  A(4  (57) 

Subtracting  the  BFP  or  (57)  from  the  total  phase  yields 

Hotal  ~  Hota,  =  A(t> flat_earth  +  A<f>z  +  A^'  “  flat _ earth  ~  AS  =  A(/>'z  (58) 
(58)  represents  the  BFPR  and  demonstrates  that  subtracting  the  BFP  from  the  total  phase 
results  in  only  the  terrain  perturbations  or  terrain  that  deviates  from  the  mean  slope 
(Figure  43). 
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The  BFPR  is  in  phase  space  and  needs  to  be  converted  to  elevation  to  be  of  any 
use  in  snow  depth  estimation.  Converting  it  to  elevation  requires  determination  of  the 
phase  to  elevation  relationship.  Recall  that  highly  accurate  elevations  were  recorded 
using  a  survey  grade  GPS  at  the  corner  reflector  locations  for  the  Mammoth  Mountain 
site.  If  the  difference  in  the  elevations  at  those  locations  to  the  average  slope  can  be 
determined,  a  direct  comparison  can  be  made  to  the  phase  in  the  BFPR  image  at  those 
same  locations.  Calculating  a  best  fit  plane  (BFP)  of  a  lower  resolution  DEM  provides 
that  average  terrain  slope.  After  the  differences  between  the  GPS  elevations  and  BFP  of 
the  DEM  are  determined,  a  linear  regression  can  be  perfonned,  which  results  in  the  phase 
to  elevation  relationship. 

The  next  step  is  to  apply  the  above  relationship  to  the  BFPR  image.  This  results 
in  an  image  that  represents  the  perturbation  of  the  terrain  from  the  average  slope.  This 
will  be  referred  to  as  the  “Best  Fit  Plane  Removed  and  Linearized”  (BFPRL)  image. 
Adding  the  BFPRL  image  to  the  previously  calculated  mean  slope  from  the  low 
resolution  DEM  results  in  a  DEM  at  the  same  resolution  as  the  SAR  image. 

The  last  step  is  to  subtract  a  Snow  Off  DEM  from  a  Snow  On  DEM  to  obtain  the 
snow  depth  on  a  per-pixel-basis. 
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or  AcJj/ 

Figure  43.  The  Best  Fit  Plane  Removed  (BFPR)  method  subtracts  the  best  fit  plane  (BFP) 
from  the  total  unwrapped  phase.  The  BFPR  isolates  the  portion  of  the  total 
unwrapped  phase  that  is  due  to  the  deviation  of  the  terrain  from  the  average  slope 
of  the  terrain  and  is  signified  as  “SAtj)”.  This  can  then  be  compared  to  ground  data 
to  establish  a  relationship  between  this  phase  and  elevation. 

1.  Best  Fit  Plane  Generation 

As  mentioned  previously,  this  method  is  applied  to  an  unwrapped  interferogram. 
From  a  data  processing  and  methodology  stand  point,  the  steps  that  were  used  in  the 
Baseline  method  starting  at  “data  preparation”  through  “phase  unwrapping”  are  the  same 
and  will  not  be  reiterated  here  and  further  discussion  will  simple  carry  on  from  that  point. 

Once  the  unwrapped  interferogram  is  available,  the  BFP  can  be  generated  from 
the  unwrapped  phase.  This  was  done  using  MATLAB.  Standard  code  to  accomplish 
orthogonal  regression  fitting  available  through  Mathworks  2013  was  used  to  generate  the 
best  fit  plane. 

Before  the  script  could  be  used,  the  2052x3646  phase  array  had  to  be  converted  to 
a  three  column  matrix.  Once  that  was  done,  the  data  were  fed  into  the  regression  script. 
The  script  itself  uses  principal  component  analysis  to  fit  a  linear  regression  that 
minimizes  the  perpendicular  distance  of  a  data  set  to  a  fitted  model,  which  is  the  BFP. 
The  output  is  in  the  fonn  of  two  vectors.  One  consists  of  two  principle  components  that 
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form  the  basis  of  a  plane.  The  other  vector  consists  of  a  principle  component  that  is 
orthogonal  to  the  first  two  components  and  defines  the  normal  to  the  plane.  In  simple 
terms,  the  code  generates  the  components  of  the  equation  of  a  plane  that  represents  the 
average  of  the  unwrapped  total  phase  by  minimizing  the  distance  of  that  plane  to  the 
phase  field. 

With  these  outputs,  the  equation  of  the  BFP  can  be  defined.  After  defining  the 
equation  it  can  be  used  to  generate  an  array  that  represents  all  the  values  on  that  plane  at 
every  pixel  location.  Once  that  array  is  established  it  can  be  used  in  subsequent  steps. 
Figure  44  is  an  example  of  the  best  fit  plane  from  the  Snow  Off  unwrapped  interferogram 
generated  from  SAR  images  21  and  22. 


Figure  44.  Best  fit  plane  generated  from  an  unwrapped  interferogram  showing  the  average 
slope  of  the  phase  image.  Mammoth  Mountain  study  site,  (37°37.7’N, 

1 19°02.7’N). 

2.  Best  Fit  Plane  Removal 

Being  able  to  generate  a  best  fit  plane  from  the  unwrapped  interferogram  is 
important  because  it  facilitates  the  isolation  of  the  phase  produced  by  high  frequency 
terrain.  As  mentioned  earlier,  this  is  accomplished  by  subtracting  the  BFP  from  the  total 
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phase.  This  effectively  removes  the  flat  earth  phase  and  the  average  slope  phase. 
Looking  back  at  (58)  this  can  be  seen. 

Hotal  -  M>total  =  flat_ earth  +  A<f>;  +  A £  ~  A<t> flat _earth  ~  =  A€  (58) 

The  actual  subtraction  of  the  best  fit  plane  is  fairly  straightforward  and  can  be 
seen  in  Figures  45  through  48.  Figure  45  is  a  representation  of  the  Snow  On  unwrapped 
interferogram  or  the  total  unwrapped  phase.  Figure  46  is  a  representation  of  the  best  fit 
plane  generated  from  that  interferogram  and  Figure  47  is  the  result  of  subtracting  46  from 
45. 


Snow  On  0MB  Urwuppod  MiOMunMnc  PMn 


Figure  45.  Unwrapped  interferogram.  Mammoth  Mountain  study  site,  (37°37.7’N, 

1 19°02.7’N). 
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Figure  46.  Best  Fit  Plane  generated  from  the  unwrapped  interferogram.  Mammoth  Mountain 

study  site,  (37°37.7’N,  1 19°02.7’N). 
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Figure  47.  This  figure  represents  the  phase  that  is  left  after  subtracting  the  BFP  from  the 
unwrapped  interferogram  and  is  referred  to  as  Best  Fit  Plane  Removed  (BFPR). 
Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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Figure  48.  Path  to  the  BFPR.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


83 


Radians 


3.  Mean  Slope  Generation  from  Low  Resolution  DEM 

The  next  step  after  the  BFPR  image  (in  phase  space)  was  established  was  to 
determine  the  relationship  between  the  phase  and  elevation.  The  first  step  in  that  process 
was  to  generate  the  elevation  slope.  This  was  accomplished  by  calculating  the  BFP  of  a 
relatively  low  resolution  DEM  of  the  same  image  area.  To  meet  this  requirement  a  10m 
resolution  DEM  was  downloaded  from  the  United  States  Department  of  Agriculture 
(USDA)  (National  Agriculture  Imagery  Program  (NAIP)  2013).  The  low  resolution 
image  was  geocorrected  using  the  four  comer  reflector  locations  and  using  a  rotation, 
scaling  translation  (RST)  transformation  to  WGS84  at  0.0858m  resolution  using  nearest 
neighbor  resampling.  Then,  one  of  the  SAR  Snow  Off  images  was  geocorrected  using 
the  corner  reflector  locations  as  previously  described  and  the  DEM  was  subset  to  match 
the  SAR  data  extent.  The  DEM  was  then  transformed  into  the  SAR  geometry  and 
converted  to  a  flat  binary  file  for  ingest  into  MATLAB  for  the  best  fit  plane  processing. 
The  resulting  DEM  can  be  seen  in  Figure  49. 


Figure  49.  3D  representation  of  a  10  meter  resolution  DEM  used  to  determine  average  slope 
of  the  terrain.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 

The  BFP,  or  average  terrain  slope,  of  the  low  resolution  DEM  can  now  be 
calculated.  The  same  process  that  was  used  to  calculate  the  BFP  of  the  unwrapped 
interferogram  was  used  in  this  case  also.  The  resulting  plane  can  be  seen  in  Figures  50. 
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Figure  50.  3D  representation  of  10m  DEM  BFP  (average  terrain  slope  in  meters).  Mammoth 

Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 

4.  Control  Point  Elevation  Difference  from  Average  Slope 

Recall  that  the  BFPR  image  from  the  unwrapped  interferogram  resulted  in  an 
image  that  represents  the  deviation  from  the  average  terrain  slope  in  phase  space.  To  get 
at  the  phase  to  height  relationship,  the  phase  perturbation  needs  to  be  compared  to  the 
elevation  perturbation.  The  GPS-derived  elevations  allow  for  a  direct  comparison 
between  the  known  elevations  at  defined  points  such  as  the  comer  reflectors  and  the 
calculated  mean  elevation  slope  (Figure  51). 

The  pixel  coordinate  for  each  GCP  was  detennined  in  the  10m  DEM  BFP  image, 
the  elevation  value  was  recorded,  and  then  compared  to  the  elevation  value  determined 
by  the  GPS.  This  difference  was  used  to  compare  the  elevation  deviation  to  the  phase 
deviation  in  the  next  step. 

5.  Phase  to  Elevation  Linear  Regression 

Richards,  2009  showed  that  the  phase  to  elevation  relationship  is  linear. 
Therefore,  direct  comparison  between  phase  slope  perturbations  and  elevation  slope 
perturbations  should  be  viable.  Figure  51  shows  the  locations  of  the  known  elevations 
with  respect  to  a  BFPR  image.  A  linear  regression  between  the  phase  and  elevation  is 
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demonstrated  in  Figure  52.  The  linear  regression  results  in  an  equation  that  can  now  be 
applied  to  the  entire  scene. 
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Figure  5 1 .  Perspective  view  of  BFPR  with  tie  points  and  accompanying  elevations  depicted. 

Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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Figure  52.  Linear  regression  for  Snow  On  01/02  SAR  image  pair  shows  the  relationship 

between  the  phase  and  elevation. 
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6.  Apply  Linear  Regression  Equation  to  BFPR  Image 

The  second  to  last  step  in  making  the  DEM  is  to  apply  the  linear  regression 
equation  to  the  BFPR  (phase)  image  seen  in  Figure  52.  The  “x”  in  the  linear  regression 
equation  represents  the  phase  array  to  which  the  equation  is  applied.  The  image  will  not 
visually  change  but  the  values  will  be  converted  from  phase  (radians)  to  elevation 
(meters).  Recall  that  this  new  image  now  represents  the  perturbation  of  the  terrain  from 
the  slope  in  elevation  or  meters. 

7.  Add  Linearized  Image  to  Average  Slope 

The  final  step  in  generating  the  DEM  is  to  add  the  above  generated  terrain 
perturbation  to  the  10m  DEM  BFP  (average  slope  (m)).  The  resulting  image  is  a  DEM  at 
the  SAR  resolution  and  can  be  seen  in  Figure  53. 


87 


0.010}  VP  UM»<< 


10m  DEM  BFP  (terrain  slope) 


Sna»On  010}  MW 


DEM 


Figure  53.  Adding  the  BFPR  after  it  has  had  the  linear  regression  equation  applied  to  the  low  resolution  DEM  BFP  (terrain  slope  in 
meters)  results  in  a  DEM  that  is  measured  in  meters.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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8.  Snow  On/Off  Differencing 

The  final  step  in  calculating  the  snow  depth  was  to  subtract  a  Snow  Off  DEM 
from  a  Snow  On  DEM.  The  resulting  image  is  the  difference  in  elevation  between  the 
two  images,  the  snow  depth,  and  can  be  seen  in  Figure  54. 
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Snow  Depth 

Figure  54.  Subtracting  the  Snow  Off  DEM  from  the  Snow  On  DEM  results  in  a  snow  depth  image.  Mammoth  Mountain  study  site, 

(37°37.7’N,  119°02.7’N). 


90 


9.  Best  Fit  Plan  Removed  Summary 

Figure  55  summarizes  the  BFPR  method. 
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Best  Fit  Plane 
Removed  Method 

Phase  (radians)  Elevation  (meters) 


Figure  55.  Best  fit  plane  removed  method  flow  chart. 
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E.  VERIFICATION 

In  order  to  determine  the  accuracy  of  the  SAR-determined  snow  depths,  manual 
snow  depth  measurements  were  taken  during  the  field  deployment.  As  mentioned  in  the 
field  work  section,  measurements  were  taken  in  a  grid  within  the  box  formed  by  the 
corner  reflectors.  Those  measurement  locations  were  recorded  with  a  GPS. 
Unfortunately  a  post  processing  survey  grade  GPS  was  not  available  at  the  time.  This 
meant  that  the  location  accuracy  was  not  the  same  as  that  of  the  corner  reflector  survey. 
The  estimated  accuracy  of  the  snow  depth  locations  was  recorded  to  be  on  the  order  of 
four  meters  in  the  horizontal. 

The  first  step  in  verifying  the  snow  depth  was  to  geo-correct  the  snow  depth 
image  using  the  previously  calculated  image-to-map  transfonnation.  Once  this  was  done, 
the  specific  locations  where  the  snow  depth  measurements  were  taken  could  be  located 
on  the  imagery  and  analyzed. 

Due  to  the  limitations  of  accuracy  of  the  GPS  locations  for  the  snow 
measurements,  an  average  snow  depth  was  calculated  for  a  radius  of  five  meters  around 
the  recorded  locations.  Figure  56  shows  a  recorded  snow  depth  image  along  with  the 
locations  of  the  manual  measurements  and  the  five  meter  radius  circles  used. 
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Figure  56.  Manual  snow  depth  measurements  were  taken  throughout  the  scene  and  are 
represented  by  the  circles  in  the  above  snow  depth  image.  The  circle  size  also 
demonstrates  the  5m  radius  area  that  was  used  to  average  the  snow  depth  in  that 
location.  Comer  reflector  locations  are  also  shown  and  labeled  (red  outlined 
circles).  Mammoth  Mountain  study  site  (312  x  176m),  (37°37.7’N,  1 19°02.7’N) 
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IV.  BASELINE  METHOD  RESULTS 


A.  PROCESSED  IMAGERY 

The  Baseline  Method  results/imagery  will  be  displayed  in  the  same  order  as 
shown  in  the  processing  flow  chart.  The  Snow  Off  case  will  be  shown  first  and  the  DEM 
produced  from  it  will  be  applied  to  the  subsequent  Snow  On  charts  for  the  calculated 
snow  depth.  Table  1  is  provided  again  for  reference. 


SAR  Image 
number 

Date 

Time 

Surface  Condition 

01 

3  April  2012 

18:47  Z,  10:47  L 

snow  covered 

02 

3  April  2012 

18:55  Z,  10:55  L 

snow  covered 

03 

3  April  2012 

19:04  Z,  11:04  L 

snow  covered 

04 

3  April  2012 

19:00  Z,  11:00  L 

snow  covered 

21 

13  July  2012 

17:58  Z,  10:58  L 

bare 

22 

13  July  2012 

18:05  Z,  11:05  L 

bare 

Table  1 .  Four  usable  Snow  On  and  two  useable  Snow  Off  SAR  images  at  0.  lm  resolution 
were  acquired  by  General  Atomics  during  the  radar  collection  phase  of  the 
research.  These  images  enabled  six  Snow  On  and  one  Snow  Off  interferometric 

pairs  to  be  calculated. 


1.  Snow  Off  Pair  21/22 

This  Snow  Off  SAR  pair  processing  started  by  first  generating  the  interferogram 
image  using  image  pairs  21  and  22.  After  the  interferogram  was  calculated  the  image 
was  then  unwrapped  (Figure  57).  The  baseline  was  then  calculated  using  both  ground 
data  and  phase  values  at  the  GCPs  (Table  2).  Both  the  range  and  azimuth  corrections 
were  calculated  and  combined  resulting  in  the  flat  earth  correction  (Figure  58).  The  flat 
earth  correction  was  subtracted  from  the  total  phase  resulting  in  the  flattened  image 
(Figures  59).  The  scaling  factor  was  then  applied  resulting  in  the  final  DEM  (Figure  60). 
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(1)  Interferogram 


Snow  Off  21/22  BL 


Snow  Off  21/22  Interferogram  Snow  Off  21/22  Unwrapped  Interferogram 


Interferogram  Unwrapped 

Interferogram 

Figure  57.  Snow  Off  SAR  image  pair  21  and  22  interferogram  and  the  resulting  unwrapped 
interferogram.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 

(2)  Calculated  Baseline.  The  system  of  equations  to  calculate 
the  baseline  between  the  two  flight  paths  requires  inputs  that  include  the  change  in  the  x, 
y,  and  z  in  meters  relative  to  the  reference  tie  point  or  comer  reflector  and  the 
interferometric  phase  at  those  tie  points.  Those  inputs  along  with  the  calculated  baseline 
can  be  seen  in  the  following  table. 
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Baseline  input  Data  21/22  (meters) 

Baseline  (meters) 

SE  (o) 

SW  (a) 

NW  (b) 

NE  (c  ) 

phaseAcf) 

-17.7481 

23.2551 

33.6268 

-3.6191 

Bx 

-9.23 

adjusted  Ax 

NA 

-92.522 

-63.6515 

28.36108 

By 

-0.7 

adjusted  Ay 

NA 

36.16798 

96.75298 

64.08585 

Bz 

123.59 

Az 

NA 

-1.813 

1.949 

4.158 

Total  BL 

123.9362 

Table  2.  The  baseline  calculation  relies  on  the  relative  A(x,  y,  z)  change  of  position  from  a 
reference  point  (SE(o))  to  three  other  known  locations  (SW(a),  NW(b),  and 
NE(c))  and  the  interferometric  phase  ( A(|) )  at  those  locations.  The  right  side  of  the 
chart  depicts  the  calculated  baseline  components  in  the  x,  y,  z  directions  (Bx,  By, 

and  Bz) 


(3)  Flat  Earth  Correction 
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Figure  58.  The  summation  of  the  range  phase  correction  and  azimuth  phase  correction  results 
in  the  total  flat  earth  phase  correction  that  will  be  applied  to  the  unwrapped 
interferogram.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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Figure  59.  Subtracting  the  flat  earth  correction  from  the  unwrapped  phase  results  in  a  phase 
that  is  entirely  due  to  the  snow  covered  terrain  and  is  known  as  the  flattened 
phase.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


(5)  Scaling  and  DEM  Generation. 

The  scaling  equation  (51)  that  converts  the  flattened  phase  to 
elevation  was  applied  to  generate  the  DEM. 
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Figure  60.  Once  the  scaling  factor  is  applied  to  the  flattened  phase,  the  resulting  image  is  a 
DEM  of  the  area  being  sampled.  The  image  on  the  right  is  the  SAR-calculated 
DEM  while  the  image  on  the  left  is  the  low  resolution  DEM  for  comparison. 

Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 

2.  Snow  On  Pair  01/02 

The  Snow  On  SAR  pair  processing  started  by  first  generating  the  interferogram 
image  using  image  pairs  01  and  02.  After  the  interferogram  was  calculated  the  image 
was  then  unwrapped  (Figure  61).  The  baseline  was  then  calculated  using  both  ground 
data  and  phase  values  at  the  GCPs  (Table  3).  Both  the  range  and  azimuth  corrections 
were  calculated  and  combined  resulting  in  the  flat  earth  correction  (Figure  62).  The  flat 
earth  correction  was  subtracted  from  the  total  phase  resulting  in  the  flattened  image 
(Figures  63).  The  scaling  factor  was  then  applied  resulting  in  the  final  DEM  (Figure  64). 
The  Snow  Off  DEM  was  then  subtracted  from  this  DEM  resulting  in  the  final  snow  depth 
image  (Figure  65). 

(1)  Interferogram 
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Figure  61.  Snow  On  SAR  image  pair  01  and  02  interferogram  and  the  resulting  unwrapped 
interferogram.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 

(2)  Calculated  Baseline.  The  system  of  equations  to  calculate 
the  baseline  between  the  two  flight  paths  requires  inputs  that  include  the  change  in  the  x, 
y,  and  z  in  meters  relative  to  the  reference  tie  point  or  comer  reflector  and  the 
interferometric  phase  at  those  tie  points.  Those  inputs  along  with  the  calculated  baseline 
can  be  seen  in  Table  3. 
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Table  3.  The  baseline  calculation  relies  on  the  relative  A(x,  y,  z)  change  of  position  from  a 
reference  point  (SE(o))  to  three  other  known  locations  (SW(a),  NW(b),  and 
NE(c))  and  the  interferometric  phase  (A(j>)  at  those  locations.  The  right  side  of  the 
chart  depicts  the  calculated  baseline  components  in  the  x,  y,  z  directions  (Bx,  By, 

and  Bz) 


(3)  Flat  Earth  Correction 
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Figure  62.  The  summation  of  the  range  phase  correction  and  azimuth  phase  correction  results 
in  the  total  flat  earth  phase  correction  that  will  be  applied  to  the  unwrapped 
interferogram.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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Figure  63.  Subtracting  the  flat  earth  correction  from  the  unwrapped  phase  results  in  a  phase 
that  is  entirely  due  to  the  snow  covered  terrain  and  is  known  as  the  flattened 
phase.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


(5)  Scaling  and  DEM  Generation. 

The  scaling  equation  was  applied  to  generate  the  DEM. 
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Figure  64.  Once  the  scaling  factor  is  applied  to  the  flattened  phase,  the  resulting  image  is  a 
DEM  of  the  area  being  sampled.  The  image  on  the  right  is  the  SAR-calculated 
DEM  while  the  image  on  the  left  is  the  low  resolution  DEM  for  comparison. 
Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(6)  Snow  Depth 
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Figure  65.  Subtracting  the  Snow  Off  DEM  from  the  Snow  On  DEM  should  result  in  an 
image  representing  the  snow  depth.  Mammoth  Mountain  study  site,  (37°37.7’N, 

1 19°02.7’N). 


3.  Snow  On  Pair  01/03 

This  Snow  On  SAR  pair  processing  started  by  first  generating  the  interferogram 
image  using  image  pairs  01  and  03.  After  the  interferogram  was  calculated  the  image 
was  then  unwrapped  (Figure  66).  The  baseline  was  then  calculated  using  both  ground 
data  and  phase  values  at  the  GCPs  (Table  4).  Both  the  range  and  azimuth  corrections 
were  calculated  and  combined  resulting  in  the  flat  earth  correction  (Figure  67).  The  flat 
earth  correction  was  subtracted  from  the  total  phase  resulting  in  the  flattened  image 
(Figures  68).  The  scaling  factor  was  then  applied  resulting  in  the  final  DEM  (Figure  69). 
The  Snow  Off  DEM  was  then  subtracted  from  this  DEM  resulting  in  the  final  snow  depth 
image  (Figure  70). 
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Figure  66.  Snow  On  SAR  image  pair  01  and  03  interferogram  and  the  resulting  unwrapped 
interferogram.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 

(2)  Calculated  Baseline.  The  system  of  equations  to  calculate 
the  baseline  between  the  two  flight  paths  requires  inputs  that  include  the  change  in  the  x, 
y,  and  z  in  meters  relative  to  the  reference  tie  point  or  comer  reflector  and  the 
interferometric  phase  at  those  tie  points.  Those  inputs  along  with  the  calculated  baseline 
can  be  seen  in  Table  4. 
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Baseline  input  Data  01/03  (meters) 

Baseline  (meters) 

SE  (o) 

SW  (a) 

NW  (b) 

NE  (c  ) 

phaseAcf) 

3.2806 

11.5804 

17.7238 

2.7155 

Bx 

0.2141 

adjusted  Ax 

NA 

-92.6249 

-63.9277 

28.17786 

By 

4.3107 

adjusted  Ay 

NA 

35.90349 

96.57073 

64.16662 

Bz 

-46.6489 

Az 

NA 

-1.813 

1.949 

4.158 

Total  BL 

46.84814 

Table  4.  The  baseline  calculation  relies  on  the  relative  A(x,  y,  z)  change  of  position  from  a 
reference  point  (SE(o))  to  three  other  known  locations  (SW(a),  NW(b),  and 
NE(c))  and  the  interferometric  phase  (A(j>)  at  those  locations.  The  right  side  of  the 
chart  depicts  the  calculated  baseline  components  in  the  x,  y,  z  directions  (Bx,  By, 

and  Bz) 


(3)  Flat  Earth  Correction 
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Figure  67.  The  summation  of  the  range  phase  correction  and  azimuth  phase  correction  results 
in  the  total  flat  earth  phase  correction  that  will  be  applied  to  the  unwrapped 
interferogram.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


(4)  Flattened  Phase  Calculation 
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Figure  68.  Subtracting  the  flat  earth  correction  from  the  unwrapped  phase  results  in  a  phase 
that  is  entirely  due  to  the  snow  covered  terrain  and  is  known  as  the  flattened 
phase.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 

(5)  Scaling  and  DEM  Generation.  The  scaling  equation  was 
applied  to  generate  the  DEM. 
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Figure  69.  Once  the  scaling  factor  is  applied  to  the  flattened  phase,  the  resulting  image  is  a 
DEM  of  the  area  being  sampled.  The  image  on  the  right  is  the  SAR-calculated 
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DEM  while  the  image  on  the  left  is  the  low  resolution  DEM  for  comparison. 
Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


(6)  Snow  Depth 
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Figure  70.  Subtracting  the  Snow  Off  DEM  from  the  Snow  On  DEM  should  result  in  an 
image  representing  or  the  snow  depth.  Mammoth  Mountain  study  site, 

(37°37.7’N,  119°02.7’N). 

4.  Snow  On  Pair  01/04 

This  Snow  On  SAR  pair  processing  started  by  first  generating  the  interferogram 
image  using  image  pairs  01  and  04.  After  the  interferogram  was  calculated  the  image 
was  then  unwrapped  (Figure  71).  The  baseline  was  then  calculated  using  both  ground 
data  and  phase  values  at  the  GCPs  (Table  5).  Both  the  range  and  azimuth  corrections 
were  calculated  and  combined  resulting  in  the  flat  earth  correction  (Figure  72).  The  flat 
earth  correction  was  subtracted  from  the  total  phase  resulting  in  the  flattened  image 
(Figures  73).  The  scaling  factor  was  then  applied  resulting  in  the  final  DEM  (Figure  74). 
The  Snow  Off  DEM  was  then  subtracted  from  this  DEM  resulting  in  the  final  snow  depth 
image  (Figure  75). 
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Figure  7 1 .  Snow  On  SAR  image  pair  0 1  and  04  interferogram  and  the  resulting  unwrapped 
interferogram.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 

(2)  Calculated  Baseline.  The  system  of  equations  to  calculate 
the  baseline  between  the  two  flight  paths  requires  inputs  that  include  the  change  in  the  x, 
y,  and  z  in  meters  relative  to  the  reference  tie  point  or  comer  reflector  and  the 
interferometric  phase  at  those  tie  points.  Those  inputs  along  with  the  calculated  baseline 
can  be  seen  in  Table  5. 
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Baseline  input  Data  01/04  (meters) 

Baseline  (meters) 

SE  (o) 

SW  (a) 

NW  (b) 

NE  (c  ) 

phaseAcf> 

5.2711 

-15.1552 

-22.2503 

-4.3794 

Bx 

5.68 

adjusted  Ax 

NA 

-92.6249 

-63.9277 

28.17786 

By 

1.45 

adjusted  Ay 

NA 

35.90349 

96.57073 

64.16662 

Bz 

-95.12 

Az 

NA 

-1.813 

1.949 

4.158 

Total  BL 

95.30047 

Table  5.  The  baseline  calculation  relies  on  the  relative  A(x,  y,  z)  change  of  position  from  a 
reference  point  (SE(o))  to  three  other  known  locations  (SW(a),  NW(b),  and 
NE(c))  and  the  interferometric  phase  ( A(|) )  at  those  locations.  The  right  side  of  the 
chart  depicts  the  calculated  baseline  components  in  the  x,  y,  z  directions  (Bx,  By, 

and  Bz) 
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Figure  72.  The  summation  of  the  range  phase  correction  and  azimuth  phase  correction  results 
in  the  total  flat  earth  phase  correction  that  will  be  applied  to  the  unwrapped 
interferogram.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


(4)  Flattened  Phase  Calculation 
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Figure  73.  Subtracting  the  flat  earth  correction  from  the  unwrapped  phase  results  in  a  phase 
that  is  entirely  due  to  the  snow  covered  terrain  and  is  known  as  the  flattened 
phase.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


(5)  Scaling  and  DEM  Generation.  The  scaling  equation  was 
applied  to  generate  the  DEM. 


Figure  74.  Once  the  scaling  factor  is  applied  to  the  flattened  phase,  the  resulting  image  is  a 
DEM  of  the  area  being  sampled.  The  image  on  the  right  is  the  SAR-calculated 
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DEM  while  the  image  on  the  left  is  the  low  resolution  DEM  for  comparison. 
Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


(6)  Snow  Depth 
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Figure  75.  Subtracting  the  Snow  Off  DEM  from  the  Snow  On  DEM  results  in  an  image 
representing  or  the  snow  depth.  Mammoth  Mountain  study  site,  (37°37.7’N, 

1 19°02.7’N). 

B.  BASELINE  METHOD  ANALYSIS 


Ultimately  the  baseline  method  proved  to  be  quite  problematic.  Based  on 
comparison  of  the  images  above,  following  the  prescribed  steps  did  not  result  in 
representative  DEMs  that  could  be  used  for  either  the  Snow  Off  or  any  of  the  Snow  On 
image  pairs  (Figures  61,  65,  70,  and  75).  The  problem  is  clearly  evident  when  the  flat 
earth  correction  is  applied.  The  resulting  image  in  phase  space  should  clearly  represent 
the  terrain  even  if  it  has  not  been  converted  from  phase  to  elevation. 

There  are  a  couple  of  potential  areas  where  the  problem  could  reside.  Confidence 
is  high  that  the  problem  does  not  reside  in  the  processing  chain  above  the  unwrapping  of 
the  interferogram.  Problems  in  those  areas  would  have  resulted  in  a  failure  to  produce 
the  fringe  images  or  a  uniform  unwrapped  interferogram.  With  the  problem  first  being 
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evident  in  the  flat  earth  removal,  the  attention  was  focused  there.  Flat  earth  removal  uses 
well  established  interferometric  theory.  It  is  therefore  unlikely  that  the  problem  is  a 
result  of  any  issues  involving  the  fundamental  equations  being  used.  Therefore  the 
problem  is  most  likely  in  the  calculations  or  the  data  that  were  used  in  the  calculations. 
The  code  that  was  used  was  scrutinized  by  both  the  author  and  outside  sources  at  Sandia 
National  labs  with  no  errors  observed.  Due  to  the  level  of  complexity  involved  in  the 
calculations,  however,  errors  in  the  coding  to  perfonn  the  calculations  cannot  be  entirely 
ruled  out.  If  it  is  assumed  that  the  calculations  were  perfonned  correctly,  the  most  likely 
source  for  the  failure  of  the  method  would  be  in  the  data  themselves. 

The  flat  earth  phase  calculation  is  composed  of  three  input  sources.  They  are  the 
incident  angle,  wavelength,  and  baseline  infonnation.  The  incident  angle  is 
straightforward  to  calculate  and  the  wavelength  is  fixed.  The  baseline  is  expected  to  be 
the  most  likely  source  of  error.  The  baseline  is  the  primary  detenninate  of  the  phase 
pattern  in  an  interferometric  image.  An  inability  to  get  an  accurate  baseline  prevents  the 
generation  of  the  flat  earth  phase  along  with  the  inability  to  derive  the  proper  phase  to 
height  ratio.  Recall  that  the  baseline  was  determined  using  a  system  of  equations  that 
used  the  phase  at  four  GCPs  and  the  relative  Ax,  y,  and  z  distances  between  those 
positions.  Errors  in  the  phase  are  not  probable  but  are  possible.  It  is  also  possible  that 
the  5  pixel  boxcar  averaging  that  was  used  for  the  corner  reflector  position  may  have 
affected  the  phase  at  that  specific  location.  This  is  not,  however,  believed  to  be  the  case. 
It  also  needs  to  be  stated  that  during  the  Snow  Off  data  collection,  the  southeast  comer 
reflector  orientation  was  affected  by  the  wind.  It  was  not  pointed  in  the  same  manner  as 
the  other  reflectors.  There  is  strong  potential  that  this  could  have  affected  the  phase  at 
that  particular  pixel  in  the  unwrapped  phase  image.  This  only  occurred  in  the  Snow  Off 
SAR  collection  and  would  not  account  for  the  inability  to  generate  reasonable  DEMs  in 
the  Snow  On  cases  as  seen  when  comparisons  were  made  to  an  independent  low 
resolution  DEM  image. 

It  is  believed  that  if  the  baseline  is  in  error,  it  is  more  likely  that  the  error  came 
from  the  A(x,y,z)  inputs  between  the  tie  points.  One  potential  source  of  error  is  the  GPS 
coordinates.  Two  measurements  were  taken  of  the  tie  point  locations  at  different  times 
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using  the  survey  grade  GPS.  After  post  processing  the  GPS  data,  the  results  were 
remarkably  consistent  with  the  locations  being  within  a  few  centimeters  from  each  other. 
This  seems  to  rule  out  tie  point  GPS  locations  as  a  source  of  error.  Another  potential 
source  of  error  with  respect  to  the  A(x,y,z)  inputs  may  involve  the  coordinate  rotation 
perfonned  to  adjust  the  real  world  A(x,y,z)  into  the  SAR  image  orientation.  The  inputs  to 
these  adjustments  were  based  on  the  aircraft’s  flight  headings  and  pointing  data  of  the 
radar.  The  limitations  of  the  accuracy  of  those  inputs  are  not  well  known  and  it  is 
difficult  to  know  how  potential  errors  in  the  recorded  flight  data  may  affect  the 
calculations  used  to  rotate  the  coordinate  system  and  thereby  cause  changes  in  the 
A(x,y,z)  inputs  for  the  baseline  calculations. 

To  explore  this  potential  problem,  sensitivity  analyses  were  perfonned  to 
detennine  the  impact  of  errors  of  A(x,y,z)  inputs  on  both  the  calculation  of  the  baseline 
components  and  the  impact  of  baseline  error  on  interferometric  phase  to  height  ratio. 
Recall  that  the  calculation  of  the  baseline  has  phase  inputs  from  the  four  different  GCPs 
along  with  the  x,  y,  and  z  change  between  the  reference  GCP  and  the  three  other  GCPs. 
Analysis  shows  that  the  overall  sensitivity  depends  on  the  length  of  the  baseline  and  can 
be  seen  in  the  table  in  Appendix  A.  This  table  was  generated  using  the  interferometric 
phase  to  height  equation  (Richards,  2009),  and  reasonable  aircraft  height,  incident  angle, 
and  wavelength  numbers  for  the  flights  performed  for  this  research.  The  analysis 
indicates  that  the  shorter  the  baseline,  the  greater  the  sensitivity  of  the  phase  to  height 
relationship.  This  is  important  because  the  flight  paths  are  designed  to  keep  the  baseline 
as  short  as  possible  to  increase  coherence.  Lynx  radar  flight  paths  typically  are  designed 
to  have  a  baselines  on  the  order  of  meters  (D.  Bickel,  2013,  personal  communication). 
Appendix  B  demonstrates  the  sensitivity  of  measurement  error  in  A(x,y,z)  measurements 
between  the  GCPs  on  the  baseline  calculation.  For  demonstration  purposes  a  baseline 
with  components  of  Bx  =  -6.14,  By  =  23.54  and  Bz  =  40.89  were  used.  This  baseline  is 
not  particularly  short  but  it  does  demonstrate  the  sensitivity  of  small  errors  in  A(x,y,z) 
inputs.  Additionally  these  charts  show  sensitivity  for  only  one  of  the  GCPs.  Remember 
each  GCP’s  A  x,  y,  and  z  affects  all  three  of  the  baseline  components,  B5jBy,  and  Bz 

differently.  Az  measurements  tend  to  be  the  most  sensitive.  Note  that  measurement 
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errors  on  the  order  of  one  centimeter  in  Az  can  result  in  a  baseline  Bz  error  of  a  quarter  of 
a  meter.  If  one  has  an  error  of  5cm  the  error  becomes  1.3  meters.  That  can  have  a 
tremendous  impact  on  the  calculated  flat  earth.  It  in  effect  adds  an  incorrect  tilt  or 
gradient  to  the  phase  field. 

The  true  sensitivity  to  potential  measurement  error  is  difficult  to  quantify  with  13 
ground  measurements  (9  A  x,  y,  z's  and  4  Acf/s)  determining  3  different  components  of  the 
baseline  (Bx  By  ,  and  Bz)  differently.  Potential  errors  ultimately  determine  the  flat  earth 
phase  that  is  subtracted  from  total  phase  and  phase  to  height  ratio.  Due  to  this 
complexity,  a  quantified  sensitivity  was  not  explored  at  this  time.  Despite  this,  it  is  clear 
from  the  examples  provided  that  these  potential  errors  make  the  task  of  determining  the 
baseline  at  the  necessary  accuracy  difficult. 

Without  the  necessary  accuracy  and  precision  in  the  ground  measurements  it  was 
impossible  to  calculate  a  baseline  that  allowed  determination  of  a  representative  flat  earth 
image  to  be  subtracted  from  the  unwrapped  interferogram.  Thus  it  was  not  possible  to 
generate  a  good  DEM  using  the  Baseline  method. 
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V.  BEST  FIT  PLANE  REMOVAL  METHOD  RESULTS 


A.  PROCESSED  IMAGERY 

The  BFPR  processed  results  will  be  broken  down  into  three  parts.  The  first  two 
parts  will  show  the  low  resolution  DEM  processing  and  the  Snow  Off  DEM  generation 
respectively.  Due  to  available  data,  only  one  Snow  Off  DEM  could  be  generated.  Both 
the  low  resolution  DEM  processing  and  Snow  Off  DEM  were  combined  with  the  Snow 
On  DEM  processed  images  to  calculate  the  snow  depth.  The  Snow  On  DEM  generation 
will  be  the  third  section  discussed  in  the  processed  imagery  portion  of  the  results. 

1.  Low  Resolution  DEM  Processing 

The  low  resolution  DEM  processing  was  straightforward.  A  10m  resolution  DEM 
was  obtained  and  geocorrected  to  match  the  observation  area  being  processed  using  the 
methods  previously  described.  Figure  76  is  a  3D  perspective  representation  of  this  DEM. 
Figure  77  is  a  3D  representation  of  the  Best  Fit  Plane  to  the  DEM.  This  BFP  was  used  in 
the  processing  of  all  of  the  higher  resolution  DEMs  for  both  the  Snow  On  and  Snow  Off 
cases. 


Figure  76.  A  relatively  low  resolution  10m  DEM  of  the  survey  area  was  obtained  from  the 
USDA  NAIP  and  then  geocorrected  to  match  the  SAR  image  area.  Centered  on 
Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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Figure  77.  Best  fit  Plane  of  low  resolution  10m  DEM  of  survey  area.  This  represents  the 
generalized  terrain  slope  for  the  Elysian  Fields  site.  Centered  on  Mammoth 
Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 

2.  Snow  Off  DEM  Generation 

(1)  Best  Fit  Plane  Generation 

The  BFP  was  generated  for  the  SAR  image  pair  21/22  by  fitting  a 
plane  to  the  total  unwrapped  phase  image  and  then  removing  it  via  subtraction  as 
previously  described  (Figure  78).  Linear  regression  was  then  used  to  convert  the  BFPR 
phase  to  the  difference  between  the  mean  slope  elevation  and  comer  reflector  location 
elevations  (Table  6,  Figure  79).  Finally,  the  10m  DEM  BFP  or  mean  elevation  slope  was 
added  back  in  to  get  the  DEM  (Figure  80). 
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Figure  78.  A  Best  Fit  Plane  (BFP)  was  generated  from  the  interferometric  SAR  image  pair  (total  unwrapped  phase)  2 1  and  22.  After 
generating  the  BFP,  it  was  subtracted  from  the  total  phase  resulting  in  a  Best  Fit  Plane  Removed  (BFPR)  image.  Mammoth 

Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(2)  Phase  Statistics 


BFPR  Statistics  Snow  Off  21/22 

X 

y 

BFP  (10m  DEM) 
(meters) 

Real  (GPS) 
(meters) 

Phase 

Difference 

(meters) 

SW  (a) 

432 

1696 

3055.95 

3064.87 

2.653742 

8.92 

NW  (b) 

678 

985 

3064.44 

3068.64 

D. 662982 

4.2 

NE  (c) 

1790 

1227 

3066.45 

3070.85 

D. 641027 

4.4 

Problem  with  SE  CR  during  collection 

Table  6.  The  phase  from  the  BFPR  image  was  compared  with  the  mean  slope  perturbation 
to  detennine  the  phase  to  elevation  relationship  for  Snow  Off  SAR  image  pairs  2 1 
and  22.  In  this  particular  case,  only  three  comer  reflectors  were  used.  The  fourth 
corner  reflector  phase  was  determined  to  be  unreliable  due  to  placement  issues 

during  the  SAR  collection 
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Figure  79.  Linear  regression  of  the  BFPR  phase  and  mean  slope  perturbation  results  in  an 
equation  that  is  applied  to  change  BFPR  image  from  phase  to  elevation  for  SAR 

image  pair  2 1  and  22. 


(3)  DEM  Generation 

At  this  point  the  linear  regression  fonnula  was  applied  to  the  BFPR 
phase  image,  resulting  in  a  BFPRL  image.  The  average  slope  was  then  added  back  in  to 
generate  the  DEM. 
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Figure  80.  Addition  of  BFPRL  to  the  10m  DEM  BFP  results  in  the  Snow  Off  SAR  image  pair  21  and  22  DEM.  Mammoth  Mountain 

study  site,  (37°37.7’N,  1 19°02.7’N). 
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Figure  8 1 .  The  BFPR  derived  DEM  (left)  generally  resembles  the  low  resolution  DEM  (right),  however,  it  exhibits  higher  detail  and  fine 

terrain  variability. 
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The  calculated  DEM  compared  well  to  the  10m  DEM  as  seen  in 
Figure  8 1 .  No  further  comparisons  of  the  calculated  DEM  to  the  low  resolution  will  be 
made  unless  there  is  a  specific  exception. 

3.  Snow  On  01/02  DEM  Generation  and  Snow  Depth  Calculation 

(1)  Best  Fit  Plane  Generation 

The  BFP  was  generated  for  the  SAR  image  pair  01/02  by  fitting  a 
plane  to  the  total  unwrapped  phase  image  and  then  removing  it  via  subtraction  as 
previously  described  (Figure  82).  Linear  regression  was  then  used  to  convert  the  BFPR 
phase  to  the  difference  between  the  mean  slope  elevation  and  corner  reflector  location 
elevations  (Table  7,  Figure  83).  Next,  the  10m  DEM  BFP  or  mean  elevation  slope  was 
added  back  in  to  get  the  Snow  On  DEM  (Figure  84).  The  last  step  was  to  subtract  the 
Snow  Off  DEM  from  the  Snow  On  DEM  to  get  the  snow  depth  (Figure  85). 
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Figure  82.  A  Best  Fit  Plane  (BFP)  was  generated  from  the  interferometric  SAR  image  pair  (total  unwrapped  phase)  01  and  02  .  After 
generating  the  BFP,  it  was  subtracted  from  the  total  phase  resulting  in  a  Best  Fit  Plane  Removed  (BFPR)  image.  Mammoth 

Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


124 


Radians 


(2)  Phase  Statistics 


BFPR  Statistics  Snow  On  01/02 

X 

y 

BFP  (10m  DEM) 
(meters) 

Real  (GPS) 
(meters) 

Phase 

Difference 

(meters) 

SW 

(a) 

432 

1696 

3055.95 

3064.87 

1.709228 

8.92 

NW 

(b) 

678 

985 

3064.44 

3068.64 

0.104561 

4.2 

NE  (c) 

1790 

1227 

3066.45 

3070.85 

-0.49067 

4.4 

SE 

1554 

1980 

3057.55 

3066.69 

1.704008 

9.14 

Table  7.  The  phase  from  the  BFPR  image  was  compared  with  the  mean  slope  perturbation 
to  determine  phase  to  elevation  relationship  for  Snow  Off  S AR  image  pairs  0 1 

and  02. 
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Figure  83.  Linear  regression  of  the  BFPR  phase  and  mean  slope  perturbation  results  in  an 
equation  that  is  applied  to  change  BFPR  image  from  phase  to  elevation  for  SAR 

image  pair  01  and02. 


(3)  DEM  Generation 

At  this  point  the  linear  regression  fonnula  was  applied  to  the  BFPR 
phase  image,  resulting  in  a  BFPRL  image.  The  average  slope  was  then  added  back  in  to 
generate  the  DEM. 
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Figure  84.  Addition  of  BFPRL  to  the  10m  DEM  BFP  results  in  the  Snow  On  SAR  image  pair  01  and  02  DEM.  The  snow  covered  DEM 

assumes  minimal  SAR  penetration.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(4)  Snow  Depth  Calculation 

The  last  step  in  the  process  was  to  subtract  the  Snow  Off  DEM 
from  the  Snow  On  DEM  to  get  the  snow  depth. 
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Figure  85.  Subtracting  the  Snow  Off  DEM  from  the  Snow  On  DEM  generated  from  the  SAR  image  pair  01  and  02  results  in  a  snow  depth 

image.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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4.  Snow  On  01/03  DEM  Generation  and  Snow  Depth  Calculation 

(1)  Best  Fit  Plane  Generation 

The  BFP  was  generated  for  the  SAR  image  pair  01/03  by  fitting  a 
plane  to  the  total  unwrapped  phase  image  and  then  removing  it  via  subtraction  as 
previously  described  (Figure  86).  Linear  regression  was  then  used  to  convert  the  BFPR 
phase  to  the  difference  between  the  mean  slope  elevation  and  comer  reflector  location 
elevations  (Table  8,  Figure  87).  Next,  the  10m  DEM  BFP  or  mean  elevation  slope  was 
added  back  in  to  get  the  Snow  On  DEM  (Figure  88).  The  last  step  was  to  subtract  the 
Snow  Off  DEM  from  the  Snow  On  DEM  to  get  the  snow  depth  (Figure  89). 
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Figure  86.  A  Best  Fit  Plane  (BFP)  was  generated  from  the  interferometric  SAR  image  pair  (total  unwrapped  phase)  1 1  and  03  .  After 
generating  the  BFP,  it  was  subtracted  from  the  total  phase  resulting  in  a  Best  Fit  Plane  Removed  (BFPR)  image.  Mammoth 

Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(2)  Phase  Statistics 


BFPR  Statistics  Snow  On  01/03 

X 

y 

BFP  (10m  DEM) 
(meters) 

R,eal(GPf  Phase 
(meters) 

Difference 

(meters) 

SW 

(a) 

432 

1696 

3055.95 

3064.87 

1.594771 

8.92 

NW 

(b) 

678 

985 

3064.44 

3068.64 

0.035803 

4.2 

NE  (c) 

1790 

1227 

3066.45 

3070.85 

-0.62206 

4.4 

SE 

1554 

1980 

3057.55 

3066.69 

1.589361 

9.14 

Table  8.  The  phase  from  the  BFPR  image  was  compared  with  the  mean  slope  perturbation 
to  determine  phase  to  elevation  relationship  for  Snow  Off  S AR  image  pairs  0 1 

and  03. 
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Figure  87.  Linear  regression  of  the  BFPR  phase  and  mean  slope  perturbation  results  in  an 
equation  that  is  applied  to  change  BFPR  image  from  phase  to  elevation  for  SAR 

image  pair  01  and  03. 


(3)  DEM  Generation 

At  this  point  the  linear  regression  fonnula  was  applied  to  the  BFPR 
phase  image,  resulting  in  a  BFPRL  image.  The  average  slope  was  then  added  back  in  to 
generate  the  DEM. 
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Figure  88.  Addition  of  BFPRL  to  the  10m  DEM  BFP  results  in  the  Snow  On  SAR  image  pair  01  and  03  DEM.  The  snow  covered  DEM 

assumes  minimal  SAR  penetration.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(4)  Snow  Depth  Calculation 

The  last  step  in  the  process  was  to  subtract  the  Snow  Off  DEM 
from  the  Snow  On  DEM  to  get  the  snow  depth. 
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Figure  89.  Subtracting  the  Snow  Off  DEM  from  the  Snow  On  DEM  generated  from  the  SAR  image  pair  01  and  03  results  in  a  snow  depth 

image.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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5.  Snow  On  01/04  DEM  Generation  and  Snow  Depth  Calculation 

(1)  Best  Fit  Plan  Generation 

The  BFP  was  generated  for  the  SAR  image  pair  01/04  by  fitting  a 
plane  to  the  total  unwrapped  phase  image  and  then  removing  it  via  subtraction  as 
previously  described  (Figure  90).  Linear  regression  was  then  used  to  convert  the  BFPR 
phase  to  the  difference  between  the  mean  slope  elevation  and  corner  reflector  location 
elevations  (Table  9,  Figure  91).  Next,  the  10m  DEM  BFP  or  mean  elevation  slope  was 
added  back  in  to  get  the  Snow  On  DEM  (Figure  92).  The  last  step  was  to  subtract  the 
Snow  Off  DEM  from  the  Snow  On  DEM  to  get  the  snow  depth  (Figure  93). 
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Figure  90.  A  Best  Fit  Plane  (BFP)  was  generated  from  the  interferometric  SAR  image  pair  (total  unwrapped  phase)  01  and  04.  After 
generating  the  BFP,  it  was  subtracted  from  the  total  phase  resulting  in  a  Best  Fit  Plane  Removed  (BFPR)  image.  Mammoth 

Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(2)  Phase  Statistics 


BFPR  Statistics  Snow  On  01/04 

X 

y 

BFP  (10m  DEM) 
(meters) 

Real  (GPS) 
(meters) 

Phase 

Difference 

(meters) 

SW 

(a) 

432 

1696 

3055.95 

3064.87 

-2.56739 

8.92 

NW 

(b) 

678 

985 

3064.44 

3068.64 

-0.5732 

4.2 

NE  (c) 

1790 

1227 

3066.45 

3070.85 

-0.78525 

4.4 

SE 

1554 

1980 

3057.55 

3066.69 

-1.01204 

9.14 

Table  9.  The  phase  from  the  BFPR  image  was  compared  with  the  mean  slope  perturbation 
to  determine  phase  to  elevation  relationship  for  Snow  Off  S AR  image  pairs  0 1 

and  04. 
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Figure  9 1 .  Linear  regression  of  the  BFPR  phase  and  mean  slope  perturbation  results  in  an 
equation  that  is  applied  to  change  BFPR  image  from  phase  to  elevation  for  SAR 

image  pair  01  and  04. 

(3)  DEM  Generation 

At  this  point  the  linear  regression  fonnula  was  applied  to  the  BFPR 
phase  image,  resulting  in  a  BFPRL  image.  The  average  slope  was  then  added  back  in  to 
generate  the  DEM. 
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Figure  92.  Addition  of  BFPRL  to  the  10m  DEM  BFP  results  in  the  Snow  On  SAR  image  pair  01  and  04  DEM.  The  snow  covered  DEM 

assumes  minimal  SAR  penetration.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(4)  Snow  Depth  Calculation 

The  last  step  in  the  process  was  to  subtract  the  Snow  Off  DEM 
from  the  Snow  On  DEM  to  get  the  snow  depth. 
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Figure  93.  Subtracting  the  Snow  Off  DEM  from  the  Snow  On  DEM  generated  from  the  SAR  image  pair  01  and  04  results  in  a  snow  depth 

image.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 


140 


Meters 


6.  Snow  On  02/03  DEM  Generation  and  Snow  Depth  Calculation 

(1)  Best  Fit  Plane  Generation 

The  BFP  was  generated  for  the  SAR  image  pair  02/03  by  fitting  a 
plane  to  the  total  unwrapped  phase  image  and  then  removing  it  via  subtraction  as 
previously  described  (Figure  94).  Linear  regression  was  then  used  to  convert  the  BFPR 
phase  to  the  difference  between  the  mean  slope  elevation  and  corner  reflector  location 
elevations  (Table  10,  Figure  95).  Next,  the  10m  DEM  BFP  or  mean  elevation  slope  was 
added  back  in  to  get  the  Snow  On  DEM  (Figure  96).  The  last  step  was  to  subtract  the 
Snow  Off  DEM  from  the  Snow  On  DEM  to  get  the  snow  depth  (Figure  98). 
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Figure  94.  A  Best  Fit  Plane  (BFP)  was  generated  from  the  interferometric  SAR  image  pair  (total  unwrapped  phase)  02  and  03.  After 
generating  the  BFP,  it  was  subtracted  from  the  total  phase  resulting  in  a  Best  Fit  Plane  Removed  (BFPR)  image.  Mammoth 

Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(2)  Phase  Statistics 


BFPR  Statistics  Snow  On  02/03 

X 

y 

BFP  (10m  DEM) 
(meters) 

R“"GPf  Phase 
(meters) 

Difference 

(meters) 

SW 

(a) 

432 

1696 

3055.95 

3064.87 

0.15168 

8.92 

NW 

(b) 

678 

985 

3064.44 

3068.64 

0.04424 

4.2 

NE  (c) 

1790 

1227 

3066.45 

3070.85 

0.09059 

4.4 

SE 

1554 

1980 

3057.55 

3066.69 

0.00938 

9.14 

Table  10.  The  phase  from  the  BFPR  image  was  compared  with  the  mean  slope  perturbation 
to  determine  phase  to  elevation  relationship  for  Snow  Off  SAR  image  pairs  02 

and  03. 
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Figure  95.  Linear  regression  of  the  BFPR  phase  and  mean  slope  perturbation  results  in  an 
equation  that  is  applied  to  change  BFPR  image  from  phase  to  elevation  for  SAR 
image  pair  02  and  03.  This  SAR  pair  has  a  particularly  low  R  value  of  .0102 
resulting  in  a  poorly  defined  DEM. 
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(3).  DEM  Generation 

At  this  point  the  linear  regression  formula  was  applied  to  the  BFPR 
phase  image,  resulting  in  a  BFPRL  image.  The  average  slope  was  then  added  back  in  to 
generate  the  DEM. 
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Figure  96.  Addition  of  BFPRL  to  the  10m  DEM  BFP  results  in  the  Snow  On  SAR  image  pair  02  and  03  DEM.  The  snow  covered  DEM 

assumes  minimal  SAR  penetration.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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In  this  particular  case  the  generated  DEM  did  not  compare  well 
with  the  low  resolution  DEM.  The  linear  regression  showed  a  very  poor  correlation 
between  the  phase  perturbation  and  the  elevation  perturbation.  The  BFPR  showed  very 
little  terrain  definition  in  the  phase  field.  This  is  most  likely  a  result  of  the  baseline 
between  this  SAR  image  pair. 

(4)  Snow  Depth  Calculation 

The  last  step  in  the  process  was  to  subtract  the  Snow  Off  DEM 
from  the  Snow  On  DEM  to  get  the  snow  depth. 


147 


M  eters 


Snow  On  DEM  Snow  Off  DEM 


Snow  Depth 

Figure  98.  Subtracting  the  Snow  Off  DEM  from  the  Snow  On  DEM  generated  from  the  SAR  image  pair  02  and  03  results  in  a  snow  depth 

image.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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7.  Snow  On  02/04  DEM  Generation  and  Snow  Depth  Calculation 

(1)  Best  Fit  Plane  Generation 

The  BFP  was  generated  for  the  SAR  image  pair  01/02  by  fitting  a 
plane  to  the  total  unwrapped  phase  image  and  then  removing  it  via  subtraction  as 
previously  described  (Figure  99).  Linear  regression  was  then  used  to  convert  the  BFPR 
phase  to  the  difference  between  the  mean  slope  elevation  and  corner  reflector  location 
elevations  (Table  11,  Figure  100).  Next,  the  10m  DEM  BFP  or  mean  elevation  slope  was 
added  back  in  to  get  the  Snow  On  DEM  (Figure  101).  The  last  step  was  to  subtract  the 
Snow  Off  DEM  from  the  Snow  On  DEM  to  get  the  snow  depth  (Figure  102). 
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Best  Fit  Plane  Removed 


Figure  99.  A  Best  Fit  Plane  (BFP)  was  generated  from  the  interferometric  SAR  image  pair  (total  unwrapped  phase)  02  and  04.  After 
generating  the  BFP,  it  was  subtracted  from  the  total  phase  resulting  in  a  Best  Fit  Plane  Removed  (BFPR)  image.  Mammoth 

Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(2)  Phase  Statistics 


BFPR  Statistics  Snow  On  02/04 

X 

y 

BFP  (10m  DEM) 
(meters) 

R“"GPf  Phase 
(meters) 

Difference 

(meters) 

SW 

(a) 

432 

1696 

3055.95 

3064.87 

4.21036 

8.92 

NW 

(b) 

678 

985 

3064.44 

3068.64 

0.59179 

4.2 

NE  (c) 

1790 

1227 

3066.45 

3070.85 

0.09139 

4.4 

SE 

1554 

1980 

3057.55 

3066.69 

2.43061 

9.14 

Table  1 1 .  The  phase  from  the  BFPR  image  was  compared  with  the  mean  slope  perturbation 
to  determine  phase  to  elevation  relationship  for  Snow  Off  SAR  image  pairs  02 

and  04. 
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Figure  100.  Linear  regression  of  the  BFPR  phase  and  mean  slope  perturbation  results  in  an 
equation  that  is  applied  to  change  BFPR  image  from  phase  to  elevation  for  SAR 

image  pair  02  and  04. 


(3)  DEM  Generation 

At  this  point  the  linear  regression  fonnula  was  applied  to  the  BFPR 
phase  image,  resulting  in  a  BFPRL  image.  The  average  slope  was  then  added  back  in  to 
generate  the  DEM. 
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Figure  101.  Addition  of  BFPRL  to  the  10m  DEM  BFP  results  in  the  Snow  On  SAR  image  pair  02  and  04  DEM.  The  snow  covered  DEM 

assumes  minimal  SAR  penetration.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(4)  Snow  Depth  Calculation 

The  last  step  in  the  process  was  to  subtract  the  Snow  Off  DEM 
from  the  Snow  On  DEM  to  get  the  snow  depth. 
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Figure  102.  Subtracting  the  Snow  Off  DEM  from  the  Snow  On  DEM  generated  from  the  SAR  image  pair  02  and  04  results  in  a  snow  depth 

image.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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8.  Snow  On  03/04  DEM  Generation  and  Snow  Depth  Calculation 

(1)  Best  Fit  Plane  Generation 

The  BFP  was  generated  for  the  SAR  image  pair  03/04  by  fitting  a 
plane  to  the  total  unwrapped  phase  image  and  then  removing  it  via  subtraction  as 
previously  described  (Figure  103).  Linear  regression  was  then  used  to  convert  the  BFPR 
phase  to  the  difference  between  the  mean  slope  elevation  and  corner  reflector  location 
elevations  (Table  12,  Figure  104).  Next,  the  10m  DEM  BFP  or  mean  elevation  slope  was 
added  back  in  to  get  the  Snow  On  DEM  (Figure  105).  The  last  step  was  to  subtract  the 
Snow  Off  DEM  from  the  Snow  On  DEM  to  get  the  snow  depth  (Figure  106). 
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Figure  103.  A  Best  Fit  Plane  (BFP)  was  generated  from  the  interferometric  SAR  image  pair  (total  unwrapped  phase)  03  and  04.  After 
generating  the  BFP,  it  was  subtracted  from  the  total  phase  resulting  in  a  Best  Fit  Plane  Removed  (BFPR)  image.  Mammoth 

Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(2)  Phase  Statistics 


BFPR  Statistics  Snow  On  03/04 

X 

y 

BFP  (10m  DEM) 
(meters) 

Real  (GPS) 
(meters) 

Phase 

Difference 

(meters) 

SW 

(a) 

432 

1696 

3055.95 

3064.87 

-4.0578 

8.92 

NW 

(b) 

678 

985 

3064.44 

3068.64 

-0.53533 

4.2 

NE  (c) 

1790 

1227 

3066.45 

3070.85 

0.552876 

4.4 

SE 

1554 

1980 

3057.55 

3066.69 

-2.09339 

9.14 

Table  12.  The  phase  from  the  BFPR  image  was  compared  with  the  mean  slope  perturbation 
to  detennine  phase  to  elevation  relationship  for  Snow  Off  SAR  image  pairs  03 

and  04. 
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Figure  104.  Linear  regression  of  the  BFPR  phase  and  mean  slope  perturbation  results  in  an 
equation  that  is  applied  to  change  BFPR  image  from  phase  to  elevation  for  SAR 

image  pair  03  and  04. 


(3)  DEM  Generation 

At  this  point  the  linear  regression  fonnula  was  applied  to  the  BFPR 
phase  image,  resulting  in  a  BFPRL  image.  The  average  slope  was  then  added  back  in  to 
generate  the  DEM. 
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Figure  105.  Addition  of  BFPRL  to  the  10m  DEM  BFP  results  in  the  Snow  On  SAR  image  pair  03  and  04  DEM.  The  snow  covered  DEM 

assumes  minimal  SAR  penetration.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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(4)  Snow  Depth  Calculation 

The  last  step  in  the  process  was  to  subtract  the  Snow  Off  DEM 
from  the  Snow  On  DEM  to  get  the  snow  depth. 
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Figure  106.  Subtracting  the  Snow  Off  DEM  from  the  Snow  On  DEM  generated  from  the  SAR  image  pair  03  and  04  results  in  a  snow  depth 

image.  Mammoth  Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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BEST  FIT  PLANE  REMOVED  ANALYSIS 


The  Best  Fit  Plane  Removed  method  showed  substantially  better  results  than  the 
baseline  method  with  respect  to  generating  viable  DEMs.  A  total  of  six  SAR  image  pairs 
were  calculated  with  varying  results  and  can  be  seen  in  Table  13.  The  snow  depth 
calculation  for  every  manual  measurement  location  and  its  associated  error  can  be  seen  in 
Appendix  C.  The  errors  were  calculated  using  two  different  data  sets.  The  first  set 
compared  the  calculated  positions  against  the  16  manually  detennined  snow  depths  in  the 
area  roughly  contained  by  the  corner  reflectors.  The  second  included  the  16  mentioned 
spots  along  with  the  five  other  measured  locations  to  the  south  of  the  box  formed  by  the 
GCPs.  The  other  five  locations  fell  within  the  edge  of  the  tree  line.  A  percentage  of 
those  calculated  snow  depths  had  radar  shadows  from  trees  in  some  of  their  averaged 
areas.  For  that  reason,  those  locations  were  not  considered  in  the  current  analysis  and 
will  be  left  for  further  consideration  at  a  future  date. 


Snow  Depth  Error  by  SAR  image  pair 

SAR  Image  Pair 

01/02 

01/03 

01/04 

02/03 

02/04 

03/04 

Average  Snow  Depth  error 
(cm) 

-8.00 

95.00 

-49.06 

175.69 

-86.56 

-41.69 

Table  13.  Average  snow  depth  error  by  SAR  image  pair. 


Examination  of  the  summary  results  in  Table  13  and  Appendix  B  shows  widely 
varying  snow  depth  errors.  Table  13  averages  all  16  snow  depth  locations  for  each  of  the 
six  SAR  image  pairs.  It  is  important  to  note  that  in  some  cases,  image  pair  01/02  in 
particular,  positive  and  negative  error  values  average  to  give  a  lower  average  error. 
Furthermore,  high  and  low  coherence  locations  are  mixed  in  these  averages.  Appendix  B 
lists  the  error  values  for  each  individual  snow  depth  location  for  each  of  the  SAR  image 
pairs. 

The  varying  snow  depth  error  results  seen  in  table  13  do  not  tell  the  whole  story. 

Each  of  the  different  SAR  image  pairs  tend  to  show  either  a  high  or  low  average  error 

rather  than  all  the  pairs  having  a  bias  in  the  same  direction.  While  it  is  not  quite  clear 

where  these  biases  are  coming  from,  it  is  believed  that  they  are  either  related  to  the  BFP 
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or  potentially  indicated  by  certain  characteristics  that  can  be  found  in  the  BFP 
calculations.  There  are  two  noted  biases.  The  first  is  the  overall  high/low  bias.  In  the 
pairs  that  have  been  computed,  two  pairs  have  a  high  bias  for  snow  depth  (01/03,  02/03), 
and  the  other  four  have  a  low  bias  (01/02,  01/04,  02/04,  03/04).  One  of  those  with  a  low 
bias,  pair  01/02,  has  only  a  slight  bias. 

Recall  that  the  flat  earth  phase  is  made  up  of  both  range  phase  and  azimuth  phase. 
It  has  also  been  observed,  at  least  in  the  cases  for  this  research,  that  the  flat  earth  phase 
tends  to  be  the  dominant  phase  in  the  unwrapped  interferograms  or  the  total  phase 
images.  This  observation  suggests  that  the  BFPs  calculated  tend  to  be  dominated  by  the 
flat  earth  phase  as  opposed  to  the  terrain.  Recall  also  that  it  is  the  flight  pair  geometry  or 
the  baseline  that  determines  the  flat  earth  phase  pattern.  A  quick  look  at  the  unwrapped 
interferograms  or  total  phase  (Figures  78,  82,  87,  92,  97,  102,  103)  show  this  to  be  true. 
Very  little  of  the  terrain  phase  can  be  seen  in  the  total  phase  images.  Often  it  is 
impossible  to  get  even  a  sense  of  the  underlying  terrain. 

Both  the  range  tilt  and  azimuth  tilt,  detennined  by  the  baseline,  can  be  observed  in 
each  BFP.  For  example,  in  SAR  Image  pair  01/02,  in  Figure  107  below,  the  range  phase 
tilts  up  toward  the  north  or  away  from  the  radar  platfonn.  In  addition  to  that  observation, 
the  azimuth  phase  also  tilts  toward  the  east  or  the  right  side  of  the  image  area.  Focusing 
on  the  azimuth  phase  and  comparing  it  to  the  high  and  low  average  errors  seem  to 
indicate  a  pattern  where  the  images  with  an  eastward  tilting  azimuth  phase  demonstrate  a 
low  average  error.  Likewise,  for  those  with  a  westward  tilting  phase  there  is  a  high 
average  error.  The  mechanism  behind  this  is  not  understood  at  this  time.  The  noted 
correlation  does  not  prove  causality;  however  it  cannot  be  ignored  and  warrants  further 
exploration. 
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01/02  8cm  Low  01/03  95cm  High  01/04  49cm  Low 


02/03  176cm  High 


02/04  87cm  Low 


03/04  42cm  Low 


Snow  On  Best  Fit  Planes 

&VP.L,  O.oi/04 


Image  pixel  position  "x"  Image  pixel  position  "x"  Image  pixel  position  "x" 


02/03 


Image  pixel  position  "x” 


02/04 


Image  pixel  position  "x" 
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Image  pixel  position  "x" 


Figure  107.  Snow  depth  error  for  each  of  the  16  manual  snow  depth  locations  for  each  of  the 
six  SAR  image  pairs  are  indicated  in  the  top  bar  graphs.  The  six  images  are  the 
BFPs  that  correspond  to  the  snow  depth  error  bar  graphs. 


The  second  bias  noted  involves  both  the  coherence  of  the  area  in  the  image  being 
looked  at  and  the  range  phase  tilt.  Note  that  there  is  a  distinct  difference  between  the  first 
eight  snow  depth  locations  and  the  second  eight  with  respect  to  coherence.  That  can 
clearly  be  seen  in  Figure  108.  The  first  8  locations  are  have  high  coherence,  ranging 
typically  above  0.85  throughout  the  SAR  image  pairs.  The  second  eight  have  lower 
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coherence,  typically  averaging  below  0.7,  with  the  majority  between  0.5  and  0.65,  and 
some  of  the  SAR  image  pairs  as  low  as  the  mid-0.3  range.  The  exact  reason  for  the  low 
coherence  in  some  of  the  SAR  image  pairs  cannot  be  absolutely  detennined.  It  is 
however  likely  that  this  is  due  to  low  values  in  the  magnitude  of  the  returns  at  those 
particular  locations.  A  quick  comparison  against  the  magnitude  images  seems  to 
corroborate  this.  Because  magnitude  is  part  of  the  coherence  calculation,  this  is  the  most 
likely  cause.  The  reason  for  the  low  return  magnitude  is  difficult  to  define  absolutely. 
The  strength  of  the  magnitude  of  the  radar  return  is  typically  due  the  surface  properties 
and  the  incidence  angle  of  the  incoming  radar  emission.  Snow  radar  reflectivity  was 
observed  to  be  significantly  less  than  that  of  the  solid  ground.  In  addition  to  that,  it 
appears  that  there  are  portions  of  the  varying  terrain  that  may  be  affected  by  a  shallow 
incidence  angle.  A  shallow  incidence  angle  could  be  responsible  for  a  decrease  in  the 
magnitude,  and  therefore  be  responsible  for  a  decrease  in  the  coherence.  One  other 
possibility  is  a  difference  in  the  liquid  water  content  of  the  snow  surface  as  opposed  to 
the  frozen  content.  Generally,  the  higher  the  liquid  content,  the  higher  the  expected 
reflectivity.  This  scenario  seems  unlikely  though  because  east  facing  slopes  should  have 
a  greater  potential  for  melting  than  west  facing  slopes  in  the  late  morning  hours  when  this 
area  was  imaged.  The  image  area  has  a  ravine  running  from  the  north  to  the  south 
through  the  scene.  The  magnitude  pattern  seen  is  the  opposite  of  what  would  be  expected 
if  there  was  disproportionate  melting  occurring.  In  other  words,  it  is  the  west  facing 
slope  that  has  the  greater  magnitudes  and  the  east  facing  slope  has  weaker  magnitudes. 
Ultimately,  the  important  take  away  is  that  there  is  a  distinct  difference  between  the  SAR 
coherency  of  the  first  and  the  second  eight  snow  depth  locations. 
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Coherence  01/02 


Figure  108.  Coherence  varies  throughout  the  scene.  The  numbered  circles  denote  the  manual 
snow  depth  validation  measurement  locations.  Note  that  measurement  locations 
1-8  have  relatively  high  coherence  compared  to  locations  9-16.  Mammoth 
Mountain  study  site,  (37°37.7’N,  1 19°02.7’N). 
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Two  approaches  were  taken  to  explore  this  potential  bias.  The  first  was  to  take  a 
look  at  the  relative  coherence  of  the  16  sites  in  all  6  SAR  image  pairs.  These  were 
plotted  along  the  “0.7”  coherence  line  that  delineated  the  two  coherence  regimes  (Figure 
108).  Note  the  pattern  that  differentiates  the  first  eight  snow  depth  locations  from  the 
second  eight  relative  to  coherence.  Normalizing  the  snow  depth  errors  for  each  of  the 
six  pairs  was  the  second  action  taken.  New  tables  were  made  with  the  average  snow 
depth  errors  subtracted  from  each  of  the  measurement  locations.  This  indicated  how  each 
individual  location  varied  with  respect  to  the  average  error.  The  results  can  be  seen  in 
Figure  109. 
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Figure  109.  Average  coherence  per  manual  snow  depth  location  for  each  of  the  six 

interferometric  image  pairs. 
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Snow  Depth  Error  (normalized  for  bias) 

01/02  01/03  01/04 


Figure  110.  This  figure  demonstrates  the  relationship  between  the  normalized  snow  depth 
error  for  each  of  the  six  interferometric  image  pairs  and  the  pattern  of  the  BFP 
derived  from  those  interferograms. 
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After  the  snow  depth  errors  were  normalized,  it  appeared  that  they  generally 
follow  the  same  pattern  as  that  seen  in  the  coherence  images  (Figure  109).  In  other 
words,  after  the  normalization,  the  first  eight  show  a  high  bias  and  the  second  show  a  low 
bias  or  the  opposite  configuration.  It  should  be  noted  that  this  signal  is  not  strong  in  all 
the  cases  such  as  SAR  image  pairs  01/04  and  03/4,  while  it  is  very  strong  in  others. 
While  the  mechanism  for  this  is  not  understood,  it  appears  that  there  may  be  some 
relationship  between  this  and  the  range  phase,  similar  to  that  which  was  seen  for  the 
azimuth  phase  relative  to  the  overall  high  or  low  average  error.  When  observing  Figure 
110,  note  that  when  the  range  phase  slopes  up  toward  the  north;  the  first  eight  snow  depth 
sites  tend  to  be  greater  than  the  nonnalized  line  while  the  second  eight  tend  to  be  less 
than  the  nonnalized  line.  When  the  range  phase  slopes  toward  the  south  or  toward  the 
radar  platform  the  opposite  is  the  case.  The  signal  does  tend  to  be  weak  in  pair  03/04  and 
nonexistent  in  image  pair  01/04.  Again  this  does  not  show  causality  but  the  consistency 
of  the  pattern  cannot  be  ignored.  Therefore  there  appears  to  be  a  potential  link  to  either 
the  coherence,  or  the  cause  of  the  coherence  pattern,  and  the  particular  flight  geometry. 

Another  observation  was  made  for  three  of  the  six  SAR  image  pairs.  Observing 
Figure  107,  note  that  in  image  pairs  01/03,  02/03,  and  02/04,  the  absolute  snow  depth 
errors  show  a  stair-step  pattern  for  every  four  measurements.  This  pattern  is  also  similar 
in  SAR  image  pair  01/02,  but  the  signal  is  not  as  strong.  Each  of  these  stairs  corresponds 
with  one  of  the  rows  in  which  the  snow  depths  were  measured  manually.  The  lower 
position  number  indicates  a  measurement  further  north  or  further  from  the  radar,  and  the 
higher  position  number  being  further  south  or  closer  to  the  antenna.  For  example,  the 
eastward  row  of  four  snow  depths  had  the  “1”  position  as  the  most  northerly  component. 
Each  successive  location  went  south  through  location  “4”  and  started  over  again  at 
position  “5”  at  the  top  of  the  collection  scene  on  the  next  row. 

In  each  of  these  four  cases  the  snow  depth  error  decreases  as  the  position  moves 
south.  This  held  true  for  every  row  regardless  of  whether  there  was  a  high  or  low  bias.  It 
also  held  true  regardless  of  the  amount  of  coherence.  There  are  a  couple  of  possibilities 
that  could  account  for  this.  The  first  one  is  that  there  may  be  an  error  in  the  overall  slope 
of  the  underlying  Snow  Off  DEM.  Recall  that  the  Snow  Off  DEM  is  subtracted  from  the 
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Snow  On  DEMs.  An  error  in  the  average  slope  of  the  Snow  Off  DEM  may  account  for 
this  pattern.  The  same  pattern  is  not  apparent  in  the  other  two  scenes,  which  raises 
questions  about  the  validity  of  this  line  of  thinking.  Another  potential  explanation  is  that 
the  error  is  contained  in  the  slope  derived  from  the  10m  DEM.  Recall  that  the  10m  DEM 
slope  was  added  back  into  both  the  Snow  On  and  Snow  Off  BFPRL  images.  If  the  slope 
has  the  wrong  tilt  it  would  be  indicated  as  an  increase  in  error  in  a  particular  direction. 
The  weakness  to  that  argument  is  that  the  same  wrong  slope  is  added  to  both  the  Snow 
On  and  Snow  Off  images.  That  should  cancel  the  error  out  when  those  images  are 
subtracted  from  each  other.  Another  potential  source  lies  with  the  BFP  generated  in  the 
Snow  On  images.  It  is  assumed  that  average  elevation  slope  for  the  Snow  On  image  is 
the  same  as  that  of  the  Snow  Off.  This  would  be  a  good  assumption  if  the  snow  laid 
evenly  across  the  scene.  We  know  that  is  not  entirely  true.  The  BFPR  images  from  the 
Snow  On  cases  may  actually  have  a  different  average  terrain  slope.  After  the  Snow  On 
BFPR  is  linearized,  it  is  added  back  in  to  the  10m  DEM  slope.  It  is  assumed  that  the 
BFPRL  image  is  a  deviation  from  the  average  slope  and  that  the  Snow  On  and  Snow  Off 
images  have  the  same  average  slope.  If  in  fact  they  don’t,  this  will  cause  a  regularly 
increasing  error  in  a  particular  direction.  For  example,  if  the  snow  depth  increases  on 
average  as  one  moves  from  the  southern  part  of  the  image  to  the  northern  part  of  the 
image,  the  snow  covered  terrain  slope  will  be  steeper  than  that  of  the  slope  calculated 
from  the  10m  DEM.  This  would  mean  that  there  would  be  in  error  in  the  slope  that  is 
added  back  in. 
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VI.  COHERENCE  OBSERVATION 


In  the  course  of  this  research  there  were  a  couple  unexpected  observations.  Recall 
that  one  of  the  reasons  that  SAR  has  not  been  used  to  determine  snow  depth  is  that  it  is 
nearly  impossible  to  get  images  with  strong  coherence  where  snow  is  involved.  This  is 
because  any  changes  in  the  surface  cause  decorrelation.  Typically,  SAR  imagery  pairs 
are  taken  over  periods  of  time  that  ensure  that  the  surface  properties  will  change.  The 
SAR  acquisitions  are  often  separated  by  at  least  several  days.  Because  of  the  design  of 
this  research  (using  multiple  acquisitions  closely  placed  in  time)  temporal  decorrelation 
was  not  an  issue.  The  acquisitions  were  taken  successively  one  after  another  insuring 
there  would  never  be  more  than  a  handful  of  minutes  between  the  acquisitions.  The 
greatest  time  difference  between  acquisitions  for  the  included  images  was  17  minutes  for 
pair  01/04.  Visually  inspecting  the  images  demonstrated  that  there  were  some 
differences  in  the  coherence.  There  is  no  discernible  pattern,  however,  until  the  images 
are  put  in  order  of  increasing  time  difference  between  the  acquisition  of  two  images  used 
in  the  generation  of  the  interferometric  pair.  This  can  be  seen  in  Figure  111. 
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INSAR  Scan  Time  Difference  and  Coherence 
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Figure  111.  Time  series  of  coherence  images  ordered  by  time  between  image  pairs.  Note  the 
decrease  in  coherency  over  time.  Mammoth  Mountain  study  site,  (37°37.7’N, 

1 19°02.7’N). 


With  the  exception  of  image  pair  03/04  with  a  time  difference  of  only  4  minutes, 
there  appears  to  be  a  very  clear  signal  of  increasing  decorrelation  between  successive 
pairs  as  the  time  difference  increases  between  the  pairs.  It  is  difficult  to  determine  the 
exact  cause  of  this  increased  decorrelation  or  why  image  pair  03/04  does  not  seem  to 
follow  the  pattern.  It  is  hypothesized  that  the  surface’s  dielectric  properties  were 
changing  rapidly  enough  that  it  could  be  seen  in  the  coherence.  Air  temperature  at  the 
collection  site  was  recorded  to  be  8°  C  at  1200  local  time  or  at  the  end  of  the  SAR 
collection  period  by  the  General  Atomics  flight  crew.  In  addition  to  that  it  was  a  sunny 
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day.  It  is  believed  that  the  surface  of  the  snow  was  increasing  its  liquid  water  content 
from  surface  melting  at  a  rate  that  was  observable  in  the  coherence  images.  It  was  also 
noted  by  the  researcher  that  at  the  beginning  of  the  SAR  data  collection  that  the  snow  was 
too  dry  to  make  a  snow  ball  that  held  together.  By  the  end  of  the  SAR  collection,  there 
was  no  difficulty  in  making  a  solid  snow  ball.  This  may  be  an  indication  that  the  liquid 
content  was  increasing  as  the  data  were  being  taken.  Unfortunately  the  time  difference 
between  snowballs  was  not  noted  and  no  snow  moisture  measurements  (or  dielectric 
measurements)  were  recorded.  It  is  believed  that  potentially  a  change  in  the  liquid  water 
content  of  the  surface  may  have  resulted  in  a  change  of  SAR  penetration  into  the  surface 
resulting  in  a  change  (decrease)  in  the  path  length,  thereby  changing  the  phase  recorded 
at  the  surface.  This  would  most  likely  result  in  decorrelation  between  images.  Greater 
differences  in  time  would  result  in  greater  changes  in  path  length  and  increased 
decorrelation. 
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VII.  CONCLUSIONS,  RECOMMENDATIONS  AND  FUTURE 

WORK 


A.  CONCLUSIONS 

The  goal  of  the  SNODAR  research  was  to  explore  the  viability  of  using  Multi¬ 
pass  Single  Look  Complex  InSAR  to  detennine  snow  depth.  Two  different  approaches 
were  taken  towards  deriving  Snow  Off  and  Snow  On  DEMs  that  could  be  differenced  to 
detennine  a  change  in  elevation  of  snow  depth.  The  first  approach  required  the 
calculation  of  the  baseline  for  the  interferometric  pair.  The  second  approach  bypassed 
the  need  for  the  baseline  knowledge  and  took  a  perturbation  or  decomposition  of  parts 
approach.  The  SAR  datasets  were  acquired  by  General  Atomics  using  a  Lynx  II  radar 
from  a  platform  flying  approximately  10km  to  the  south  of  the  study  area. 

1.  Baseline  Method 

The  Baseline  method  approach  ultimately  ended  up  being  problematic.  The 
method  attempted  to  detennine  the  baseline  between  SAR  collections  by  solving  a  system 
of  equations  using  ground  control  points  and  their  associated  phases.  The  DEMs 
generated  from  this  approach  were  ultimately  unrepresentative  of  the  real  terrain.  The 
problem  revolved  around  the  ability  to  get  a  correct  baseline.  This  was  detennined  by 
observing  that  the  removal  of  the  flat  earth  phase  from  interferogram  did  not  result  in  a 
match  to  known  terrain  characteristics.  The  flat  earth  phase  is  wholly  dependent  on  an 
accurate  knowledge  of  the  baseline.  Sensitivity  tests  indicate  that  there  was  very  little 
room  for  input  error  when  calculating  the  baseline. 

The  sensitivity  of  the  Baseline  method  to  measurement  error  is  a  challenging 
problem.  Limitations  in  the  ability  to  accurately  quantify  the  input  data  will  be  difficult 
to  overcome.  Changing  the  baseline  configuration  of  the  flight  paths  to  increase  the 
baseline  will  decrease  the  sensitivity  to  potential  measurement  error.  There  are,  however, 
limits  to  the  size  of  the  baseline  before  coherence  becomes  an  issue.  In  addition,  there 
may  be  approaches  that  will  improve  the  input  data  for  the  system  of  equations  used  to 
calculate  the  baseline.  The  baseline  calculations  proved  to  be  most  sensitive  to  the 
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vertical  component  or  Az  difference  between  the  ground  control  points.  It  was  shown 
that  small  potential  errors  could  result  in  large  errors  in  the  flat  earth  phase  calculations. 

This  method  was  mathematically  rigorous  and  complicated  both  from  a 
theoretical  point  of  view  and  a  computational  point  of  view.  This  complexity  interjected 
many  opportunities  for  either  error  in  the  processing  or  error  in  the  computations.  It  is 
not  clear  from  this  research  why  the  Baseline  method  did  not  work  in  this  case,  and  it  is 
believed  that  the  fundamentals  are  still  there  for  a  successful  application  of  the  method. 
This  is  a  topic  for  future  research — several  suggestions  are  offered  below  and  in  the 
Recommendations  section. 

2.  Best  Fit  Plane  Removed  Method 

The  Best  Fit  Plane  Removed  method  showed  greater  success  relative  to  the 
baseline  method.  The  BFPR  method  bypassed  the  requirement  for  the  baseline 
knowledge  by  taking  a  perturbation  or  decomposition  approach.  This  was  done  by 
generating  a  best  fit  plane  to  the  interferometric  image  and  then  removing  it  from  the 
total  phase  data.  This  isolated  the  interferometric  phase  caused  by  the  terrain  that 
deviated  from  the  mean  slope.  It  also  removed  the  flat  earth  phase  that  was  an  issue  with 
the  Baseline  method.  The  phase  was  compared  to  the  difference  of  the  elevation  at  GCPs 
relative  to  the  average  elevation  slope.  The  average  slope  was  determined  using  a 
relatively  low  resolution  DEM  to  generate  the  mean  slope  and  comparing  it  to  survey 
grade  GPS  elevations.  A  linear  regression  was  performed  and  applied  to  the  “Best  Fit 
Plane  Removed”  image  to  convert  phase  to  elevation.  That  image  was  then  added  back 
to  the  average  slope,  resulting  in  a  DEM.  A  Snow  Off  DEM  compared  well  to  the  real 
world  10m  DEM.  After  computing  DEMs  from  both  Snow  On  and  Snow  Off  scenes  they 
were  differenced  to  calculate  snow  depth. 

The  snow  depth  results  of  six  Snow  On  SAR  pairs  were  compared  to  16  manually 
measured  snow  depth  locations  with  varying  degrees  of  success.  The  SAR  image  pairs 
showed  an  average  error  of  -8cm,  95cm,  -49cm,  175cm  87cm  and  42cm  for  the 
respective  six  SAR  pairs.  The  results  also  indicated  that  coherence  of  the  unwrapped 
InSAR  image  played  a  role  in  the  DEM  generation.  Of  the  16  locations  manually 
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measured,  eight  fell  in  a  high  coherence  regime  indicated  by  coherences  greater  than  0.7 
and  the  others  fell  in  a  regime  indicated  by  coherence  less  than  0.7.  In  almost  all  of  the 
cases  the  magnitude  of  the  error  for  each  of  the  SAR  image  pairs  fell  in  two  categories 
determined  by  this  regime. 

There  did  appear  to  be  a  consistent  pattern  of  either  high  or  low  bias  in  the  BFPR- 
calculated  snow  depth  results.  Four  of  the  SAR  image  pairs  demonstrated  a  low  average 
for  the  snow  depths  while  the  other  two  pairs  demonstrated  a  high  average.  The  pattern 
indicates  that  errors  may  be  either  related  to  or  indicated  by  the  BFPs  produced  from  the 
unwrapped  interferograms.  There  appear  to  be  two  different  biases.  The  first  is  that  the 
slope  of  the  azimuth  aspect  of  the  BFP  affects  the  direction  of  the  bias.  It  was  observed 
that  an  eastward  tilt  in  the  BFP  was  consistent  with  SAR  pairs  with  a  bias  towards  low 
snow  depth  errors.  Those  with  a  westward  tilt  demonstrated  a  bias  towards  high  snow 
depth  errors.  The  second  bias  is  not  as  well  defined  but  does  appear  to  exist  in  this 
dataset.  This  bias  appears  to  relate  coherence  in  the  image  and  the  range  slope  of  the 
BFP.  After  nonnalizing  the  error  there  was  a  clear  difference  between  the  snow  depth 
locations  with  high  and  low  coherence.  Additionally,  the  determination  of  whether  the 
high  or  low  coherence  was  above  or  below  the  nonnalization  line  appeared  to  be 
detennined  by  the  range  tilt  of  the  BFP.  This  pattern  is  not  fully  understood. 
Furthennore,  the  observed  pattern  does  not  necessarily  indicate  causality.  Additional 
SAR  image  pairs  should  be  tested  to  ensure  that  the  pattern  holds. 

Another  observation  was  made  in  four  of  the  six  SAR  image  pairs.  It  appeared 
that  regardless  of  the  coherence,  the  calculated  error  decreased  as  the  observations  moved 
southward  or  in  the  direction  toward  the  sensor.  This  is  indicative  of  an  issue  in  the  slope 
of  one  or  more  of  the  BFPR  elements.  Slope  issues  could  arise  from  the  calculation  of 
the  BFP,  accuracy  of  the  low  resolution  DEM  used  to  detennine  the  deviation  of  the  high 
frequency  terrain  from  the  average  slope,  or  an  issue  with  representativeness  of  the  low 
resolution  DEM  relative  to  the  true  slope  of  the  snow  covered  terrain. 

Furthennore  this  research  demonstrates  that  the  Ku-band  radar  is  capable  of 

discerning  the  snow  air  interface  (minimal  penetration).  This  is  evident  in  both  its  ability 

to  see  features  on  the  snow  surface  such  as  tracks  in  the  snow  from  the  researchers,  and 
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the  representative  DEMs  that  consistently  showed  the  terrain  or  snow  surface.  The 
DEMs  also  consistently  showed  the  Snow  On  terrain  to  be  higher  than  that  of  the  Snow 
Off  terrain. 


3.  Coherence  Observation 

One  final  observation  was  unexpected.  It  was  assumed  from  the  beginning  of  this 
research  that  the  decorrelation  normally  associated  with  InSAR  pairs  from  more 
traditional  collections  with  time  separation  of  multiple  days  would  be  nonexistent.  The 
greatest  time  difference  between  SAR  image  pairs  was  17  minutes  with  the  shortest  time 
difference  being  4  minutes.  It  was  observed  that  when  the  coherence  images  were 
ordered  according  to  increasing  time  between  the  images,  that  the  coherence  decreased 
from  one  image  to  the  next.  This  indicated  that  decorrelation  was  occurring  on  a  very 
short  time  scale.  This  suggests  that  the  dielectric  properties  of  the  snow  surface  were 
rapidly  changing  during  the  SAR  data  acquisition.  This  is  only  speculative,  however, 
since  no  specific  snow-moisture  or  other  measurements  were  taken  of  the  surface  as  a 
function  of  time. 

4.  Implications 

While  perfect  results  were  not  achieved,  the  BFPR  method  shows  promise.  The 
ground  work  has  been  laid  for  further  investigation.  In  addition,  while  penetration  into 
the  snow  pack  by  the  Ku-band  radar  was  not  directly  tested,  it  appears  that  it  was 
minimal.  The  greatest  challenge  in  this  research  was  achieving  good  DEMs  utilizing 
multiple  SAR  passes  with  an  aircraft  with  only  one  antenna  and  an  unknown  baseline. 
SAR  acquisition  using  an  aircraft  equipped  with  a  bistatic  antenna  system  with  a 
frequency  in  the  Ku-band  or  higher  would  greatly  simplify  the  process  and  increase  the 
probability  of  successful  snow  depth  detennination.  While  one  of  the  main  goals  was  to 
derive  a  method  that  could  be  used  with  operational  monostatic  platforms,  it  would 
benefit  future  research  to  test  these  techniques  with  a  system  that  is  better  suited  for 
making  DEMs.  The  nature  of  the  technique  would  not  then  be  dependent  on  the  ability  to 
derive  DEMs  using  monostatic  SAR  platforms  with  the  attendant  baseline  detennination 
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problems.  Once  SAR  interaction  with  the  snow  surfaces  is  better  codified;  the  focus 
could  transition  to  the  platforms  with  only  one  antenna. 

B.  RECOMMENDATIONS  AND  FUTURE  WORK 

1.  Baseline  Recommendations 

While  the  Baseline  method  proved  to  be  problematic,  there  are  still  some  potential 
things  that  can  be  done  to  move  it  toward  success.  Changing  the  flight  parameters  such 
as  lengthening  baseline  should  reduce  sensitivity  to  any  introduced  measurement  error. 
Additionally,  the  placement  of  the  GCPs  in  a  manner  that  increases  the  difference  in  the 
vertical  should  minimize  ambiguity  caused  by  the  CR  elevations  being  too  close  together. 

2.  Explore  the  causes  of  the  Best  Fit  Plane  Removed  Biases 

The  BFPR  method  has  shown  potential  but  still  has  aspects  that  are  not  fully 
understood.  The  noted  biases  appeared  to  either  be  caused  by  or  indicated  by  the  BFP 
determined  from  the  unwrapped  interferogram.  While  the  bias  pattern  was  quite 
consistent  across  all  six  SAR  image  pairs,  it  does  not  prove  causality.  The  mechanisms 
are  not  understood.  Collection  and  analysis  of  additional  data  sets  would  help  determine 
if  this  bias  to  BFP  relationship  holds.  In  addition,  once  these  biases  are  understood,  it 
should  lead  to  increased  accuracy. 

3.  Circumvent  the  biases 

If  the  hypothesized  biases  hold  to  be  true,  there  may  be  a  way  to  circumvent  them. 
Recall  that  the  flat  earth  phase  pattern  is  determined  by  the  baseline  of  the  aircraft  pair 
flight  geometry.  The  interferograms  are  generated  from  a  master/slave  relationship.  If 
this  relationship  is  reversed,  the  resulting  interferogram  will  have  a  reversed  BFP  in  both 
the  range  and  azimuth  directions.  If  this  happens  and  the  biases  are  related  to  SAR 
geometry,  one  should  be  able  get  the  exact  opposite  bias  when  the  data  are  processed. 
The  results  from  the  opposite  master/slave  relationships  with  the  presumed  opposite 
biases  should  be  able  to  be  averaged  for  the  correct  solution. 
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4. 


Bistatic  Antenna 


The  purpose  of  this  research  was  to  determine  if  snow  depth  could  be  detennined 
from  SAR  interferometry.  Generating  a  good  DEM  was  the  greatest  challenge  in  this 
research  because  of  limitations  imposed  by  the  use  of  the  monostatic  (single  antenna) 
SAR  system.  The  monostatic  configuration  is  commonly  used  for  SAR  imagery 
acquisition  and  those  data  were  what  was  available  for  this  research.  An  aircraft  with  a 
bistatic  antenna  (dual  antenna)  configuration  would  be  much  better  suited  for  DEM 
generation.  The  baseline  would  always  be  known  to  a  high  degree  of  precision. 
Knowing  the  baseline  would  make  much  of  the  DEM  generation  significantly  easier  and 
allow  for  the  focus  of  the  research  to  be  on  the  SAR  interaction  with  the  snow. 

5.  Negation  of  ground  data 

To  be  a  viable  technique,  this  capability  will  eventually  need  to  be  able  to  be 
performed  without  highly  quantified  GCPs.  Different  methods  should  be  explored  to 
determine  if  this  is  possible. 

6.  Snow  conditions 

Snow  conditions  affect  the  ability  to  get  returns  from  the  surface  and  have  the 
potential  to  affect  the  overall  accuracy.  Varying  snow  conditions  from  different  times  of 
the  season  with  different  properties  should  be  explored  to  detennine  the  effects  on  this 
technique. 


7.  Collection  method 

There  are  several  aspects  of  the  collection  that  should  be  explored.  The  baseline 
configuration  as  mentioned  above  may  play  a  role  in  the  BFPR  method.  Additionally 
collecting  data  using  a  strip-map  approach  to  cover  large  areas  of  land  should  be 
considered  as  should  using  differing  resolutions. 

8.  Perturbation  isolation 

The  BFPR  approach  was  used  to  isolate  the  terrain  phase  perturbation  for  this 
research.  Other  approaches  are  possible.  For  example  a  fast  Fourier  transfonn  (FFT) 
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could  be  used  to  separate  the  low  frequency  terrain  slope  from  the  higher  frequency 
terrain  detail  (Shurong  et  al.  2009).  Examination  of  this  FFT  approach  may  provide 
additional  insight  to  the  BFPR  results. 
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APPENDIX  A.  BASELINE  TO  M/RAD  SENSITIVITY 


Based  on  Eq  6.10 


Sensitivity  of  M/Rad  to  change  in  Baseline  (Potential  Error  analysis) 


baseline 

in 

meters 

0 

degrees 

lambda 

height  in 
meters 

Pto  H 

M/Rad 

1 

77 

0.018 

2334.5 

0.0155 

64.3879 

2 

77 

0.018 

2334.5 

0.0311 

32.1940 

3 

77 

0.018 

2334.5 

0.0466 

21.4626 

4 

77 

0.018 

2334.5 

0.0621 

16.0970 

5 

77 

0.018 

2334.5 

0.0777 

12.8776 

6 

77 

0.018 

2334.5 

0.0932 

10.7313 

7 

77 

0.018 

2334.5 

0.1087 

9.1983 

8 

77 

0.018 

2334.5 

0.1242 

8.0485 

9 

77 

0.018 

2334.5 

0.1398 

7.1542 

10 

77 

0.018 

2334.5 

0.1553 

6.4388 

11 

77 

0.018 

2334.5 

0.1708 

5.8534 

12 

77 

0.018 

2334.5 

0.1864 

5.3657 

13 

77 

0.018 

2334.5 

0.2019 

4.9529 

14 

77 

0.018 

2334.5 

0.2174 

4.5991 

15 

77 

0.018 

2334.5 

0.2330 

4.2925 

16 

77 

0.018 

2334.5 

0.2485 

4.0242 

17 

77 

0.018 

2334.5 

0.2640 

3.7875 

18 

77 

0.018 

2334.5 

0.2796 

3.5771 

19 

77 

0.018 

2334.5 

0.2951 

3.3888 

20 

77 

0.018 

2334.5 

0.3106 

3.2194 

21 

77 

0.018 

2334.5 

0.3261 

3.0661 

22 

77 

0.018 

2334.5 

0.3417 

2.9267 

23 

77 

0.018 

2334.5 

0.3572 

2.7995 

24 

77 

0.018 

2334.5 

0.3727 

2.6828 

25 

77 

0.018 

2334.5 

0.3883 

2.5755 

26 

77 

0.018 

2334.5 

0.4038 

2.4765 

27 

77 

0.018 

2334.5 

0.4193 

2.3847 

28 

77 

0.018 

2334.5 

0.4349 

2.2996 

29 

77 

0.018 

2334.5 

0.4504 

2.2203 

30 

77 

0.018 

2334.5 

0.4659 

2.1463 

31 

77 

0.018 

2334.5 

0.4815 

2.0770 

32 

77 

0.018 

2334.5 

0.4970 

2.0121 

183 

1  meter 
change  in 
base  line 


Difference  in  Height  (z)  to  Phase 
(A<|>)  ratio,  M/Rad  change 


1 

to 

2 

32.194 

2 

to 

3 

10.731 

3 

to 

4 

5.366 

4 

to 

5 

3.219 

5 

to 

6 

2.146 

6 

to 

7 

1.533 

7 

to 

8 

1.150 

8 

to 

9 

0.894 

9 

to 

10 

0.715 

10 

to 

11 

0.585 

11 

to 

12 

0.488 

12 

to 

13 

0.413 

13 

to 

14 

0.354 

14 

to 

15 

0.307 

15 

to 

16 

0.268 

16 

to 

17 

0.237 

17 

to 

18 

0.210 

18 

to 

19 

0.188 

19 

to 

20 

0.169 

20 

to 

21 

0.153 

21 

to 

22 

0.139 

22 

to 

23 

0.127 

23 

to 

24 

0.117 

24 

to 

25 

0.107 

25 

to 

26 

0.099 

26 

to 

27 

0.092 

27 

to 

28 

0.085 

28 

to 

29 

0.079 

29 

to 

30 

0.074 

30 

to 

31 

0.069 

31 

to 

32 

0.065 

32 

to 

33 

0.061 

33 

77 

0.018 

2334.5 

0.5125 

1.9511 

34 

77 

0.018 

2334.5 

0.5280 

1.8938 

35 

77 

0.018 

2334.5 

0.5436 

1.8397 

36 

77 

0.018 

2334.5 

0.5591 

1.7886 

37 

77 

0.018 

2334.5 

0.5746 

1.7402 

38 

77 

0.018 

2334.5 

0.5902 

1.6944 

39 

77 

0.018 

2334.5 

0.6057 

1.6510 

40 

77 

0.018 

2334.5 

0.6212 

1.6097 

41 

77 

0.018 

2334.5 

0.6368 

1.5704 

42 

77 

0.018 

2334.5 

0.6523 

1.5330 

43 

77 

0.018 

2334.5 

0.6678 

1.4974 

44 

77 

0.018 

2334.5 

0.6834 

1.4634 

45 

77 

0.018 

2334.5 

0.6989 

1.4308 

46 

77 

0.018 

2334.5 

0.7144 

1.3997 

47 

77 

0.018 

2334.5 

0.7300 

1.3700 

48 

77 

0.018 

2334.5 

0.7455 

1.3414 

33 

to 

34 

0.057 

34 

to 

35 

0.054 

35 

to 

36 

0.051 

36 

to 

37 

0.048 

37 

to 

38 

0.046 

38 

to 

39 

0.043 

39 

to 

40 

0.041 

40 

to 

41 

0.039 

41 

to 

42 

0.037 

42 

to 

43 

0.036 

43 

to 

44 

0.034 

44 

to 

45 

0.033 

45 

to 

46 

0.031 

46 

to 

47 

0.030 

47 

to 

48 

0.029 

48 

to 

49 

Table  14.  The  left  table  demonstrates  the  sensitivity  of  the  height  to  phase  ratio  in  an 
interferogram  based  on  the  change  in  the  baseline  length.  The  right  table 
demonstrates  the  actual  change  height  to  phase  relationship  from  1  meter  in 

baseline  to  the  next. 
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Figure  1 12.  Baseline  length  to  interferometric  phase  factor  relationship  based  on  Table  14. 
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APPENDIX  B.  BASELINE  COMPONENT  SENSITIVITY  TO 
GROUND  DIFFERENCE  ERROR 
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Figure  113.  This  figure  demonstrates  the  change  that  the  input  values  of  the  A  x,  y,  and  z 

elements  into  the  baseline  calculation  system  of  equations  has  on  the  length  of  the 

calculated  baseline. 
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APPENDIX  C.  CALCULATED  SNOW  DEPTH 


Calcu 

lated  Snow  depth  (cm) 

Manual 

01/02 

01/03 

01/04 

02/03 

02/04 

03/04 

1 

175 

220 

339 

206 

626 

-13 

156 

2 

164 

168 

270 

125 

428 

-28 

112 

3 

164 

169 

259 

131 

349 

-6 

120 

4 

130 

144 

224 

86 

235 

-16 

112 

5 

184 

196 

361 

135 

634 

-13 

131 

6 

183 

177 

327 

97 

527 

20 

126 

7 

164 

183 

308 

119 

463 

47 

149 

8 

239 

210 

331 

115 

431 

111 

172 

9 

209 

198 

353 

215 

620 

161 

164 

10 

180 

160 

268 

195 

488 

185 

137 

11 

166 

138 

224 

170 

330 

195 

133 

12 

154 

120 

203 

98 

188 

151 

111 

13 

169 

149 

242 

88 

171 

117 

86 

14 

135 

110 

176 

69 

43 

142 

109 

15 

136 

108 

167 

18 

-19 

117 

100 

16 

130 

104 

150 

30 

-21 

127 

97 

Off  grid  snow  depths  (trees  present) 

17... 18 

126 

216 

202 

79 

78 

157 

163 

18... 19 

133 

167 

195 

96 

165 

0.89 

169 

19.. .20 

138 

115 

208 

59 

312 

-29 

163 

20.. .21 

133 

87 

152 

146 

313 

173 

134 

21.. .22 

162 

109 

149 

109 

43 

167 

139 

Table  15.  Calculated  snow  depth  for  each  location  for  each  of  the  six  SAR  image  pairs  for 

Best  Fit  Plane  Removed  Method. 
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APPENDIX  D.  SNOW  DEPTH  ERROR 


Error  in  Snow  Depth  Calculated  (cm) 

01/02 

01/03 

01/04 

02/03 

02/04 

03/04 

1 

45 

164 

31 

451 

-188 

-19 

2 

4 

106 

-39 

264 

-192 

-52 

3 

5 

95 

-33 

185 

-170 

-44 

4 

14 

94 

-44 

105 

-146 

-18 

5 

12 

177 

-49 

450 

-197 

-53 

6 

-6 

144 

-86 

344 

-163 

-57 

7 

19 

144 

-45 

299 

-117 

-15 

8 

-29 

92 

-124 

192 

-128 

-67 

9 

-11 

144 

6 

411 

-48 

-45 

10 

-20 

88 

15 

308 

5 

-43 

11 

-28 

58 

4 

164 

29 

-33 

12 

-34 

49 

-56 

34 

-3 

-43 

13 

-20 

73 

-81 

2 

-52 

-83 

14 

-25 

41 

-66 

-92 

7 

-26 

15 

-28 

31 

-118 

-155 

-19 

-36 

16 

-26 

20 

-100 

-151 

-3 

-33 

Average  Grid 

-8.00 

95.00 

-49.06 

175.69 

-86.56 

-41.69 

Off  grid  snow  depths  (trees  present) 

17... 18 

90 

76 

-47 

-48 

31 

37 

18... 19 

34 

62 

-37 

32 

-132.11 

36 

19.. .20 

-23 

70 

-79 

174 

-167 

25 

20.. .21 

-46 

19 

13 

180 

40 

1 

21. ..22 

-53 

-13 

-53 

-119 

5 

-23 

Average  all 

-6.00 

82.57 

-47.05 

144.29 

-76.58 

-28.14 

Table  16.  Calculated  error  in  the  snow  depth  for  each  location  for  each  of  the  six  SAR 

image  pairs  for  the  Best  Fit  Plane  Removed  method. 
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